key: cord-0171538-kirbycpf authors: Ruan, Guangchun; Wu, Dongqi; Zheng, Xiangtian; Sivaranjani, S.; Zhong, Haiwang; Dahleh, Munther A.; Kang, Chongqing; Xie, Le title: A Cross-Domain Approach to Analyzing the Short-Run Impact of COVID-19 on the U.S. Electricity Sector date: 2020-05-11 journal: nan DOI: nan sha: fae46bbb747ebfffa149fab21467e5797a3a9e49 doc_id: 171538 cord_uid: kirbycpf The novel coronavirus disease (COVID-19) has rapidly spread around the globe in 2020, with the U.S. becoming the epicenter of COVID-19 cases since late March. As the U.S. begins to gradually resume economic activity, it is imperative for policymakers and power system operators to take a scientific approach to understanding and predicting the impact on the electricity sector. Here, we release a first-of-its-kind cross-domain open-access data hub, integrating data from across all existing U.S. wholesale electricity markets with COVID-19 case, weather, cellular location, and satellite imaging data. Leveraging cross-domain insights from public health and mobility data, we uncover a significant reduction in electricity consumption across that is strongly correlated with the rise in the number of COVID-19 cases, degree of social distancing, and level of commercial activity. . Architecture of COVID-EMDA + data hub. (a) Processing flowchart of COVID-EMDA + data hub. Heterogeneous data sources are handled, including electricity market, COVID-19 cases, cellular location data, satellite imagery and weather data. To coordinate five data sources in the same geographical scales, the geocoding technique is applied to transform COVID-19 cases and weather data. The entire processing reflects the objective of data consistence, data compaction and data checking. (b) The architecture contains the date, data category and location dimensions. The main dimension is the date dimension due to the importance of time-series relationships. Along the main dimension, one can retrieve multiple data slices or spreedsheet data files. The yellow cubic represents one of the load datasets in New York City. (c) Map of the United States representing the regions of operation of the seven RTOs or electricity markets. A cross-market comparison, with both the point-and interval-estimation results, is conducted in Tab. 1 to show the impact of COVID-19 on different marketplaces. The interval estimation is calculated using the 10% and 90% quantiles which can be regarded as reliable estimation boundaries. The backcast models successfully capture the dynamics of changes in electricity consumption and provide a statistical comparison among different markets. It is clearly seen that all the markets experienced a reduction in electricity consumption in April; however, the magnitudes of the reductions are diverse, varying from 6.71% to 11.92%. The dynamics of changes in electricity consumption are also varied -the demand in NYISO, ISO-NE and CAISO started dropping in March, while other markets remained as yet unchanged (note that interval estimations contain zero). Additionally, our estimation results for April match well with official reports [10] [11] [12] . According to Table 1 , in April, NYISO and ISO-NE experienced the most severe reduction in electricity consumption, (11.92% and 11.28% respectively), while ERCOT and SPP suffered the least. In dense urban areas, the impact was more pronounced, with New York City and Boston experiencing a 15.37% and 13.70% reduction in electricity consumption respectively in April, likely due to the high population density and large share of commercial energy use in these areas. (The same factors explain why cities like Houston and Los Angeles that are more geographically dispersed were not significantly impacted). We will examine these factors more closely in the following section. An increase in the number of COVID-19 cases indirectly impacts electricity consumption as it results in isolation policies which in turn lead to an increase in social distancing (size of the stay-at-home population), as well as shut down of businesses (commercial loads). This trend is clearly discernible in cellular location data as an increase in stay-at-home population (Supplementary Fig. S -2) and a reduction in visits to retail establishments (Supplementary Fig. S-3 ). Fig. 4 shows the multi-dimensional relationship between the number of COVID-19 cases, social distancing, shut down rate of commercial activity, and electricity consumption. In order to observe the impact of the rise in both COVID-19 cases and social distancing and shut down of commercial activity on the reduction in electricity consumption, we apply another series of backcast models to estimate daily average electricity consumption and calculate the time-varying reduction rate. Fig. 9 shows the trace of the evolution of daily new confirmed cases and social distancing, and the associated rate of reduction in electricity consumption for two representative metropolises -New York City and Philadelphia, indicating a fast developing period in March 2020 and a more stable period afterwards. Similar dynamics are observed in other COVID-19 hotspot cities ( Supplementary Fig. S-4) . The trace of evolution of electricity consumption demonstrates the dynamically evolving, multi-dimensional relationship between the number of daily new COVID-19 cases, the size of the stay-at-home population, and the reduction in electricity consumption. A unique insight from the time-series description of such cross-domain data is the wide variation in the time scales of the changes in the electricity consumption, public health, stay-at-home, work-on-site, and retail mobility data (Fig. 6 ). The mobility in the retail sector has the earliest response in terms of the rate of change (gradually dropping from late February 2020 and continuing to go down until late April 2020), resulting from bottom-up responses of consumers to the emerging pandemic. On the other hand, the population of on-site workers shows a sharp, abrupt change right around mid-March, as a result of top-down federal and state-level policy decisions such as stay-at-home orders. The electricity consumption shows a slightly delayed reduction with respect to the number of COVID cases. This insight, to our best knowledge, is first revealed in Fig. 6 , and suggests a very different efficacy of social distancing arising from top-down government policies and from bottom-up individual responses. We now rigorously quantify the multi-dimensional relationship shown in Fig. 4 30 model. Such VAR models have been widely adopted in econometrics 31 , and have also been explored in the context of electricity markets 32 . The VAR model captures how the first order differential of the reduction in electricity consumption on day t is correlated with different factors, namely, the daily new confirmed COVID-19 cases, the stay-at-home population and the population of on-site workers (indicative of social distancing), as well as the mobility in the retail sector (indicative of commercial electricity loads), on days t − i, where i = 1, . . . , p, and p is the user-defined order of the VAR model (see the Methods section for the definition, and Supplementary Methods SM-1,2,3 for details of the calibration, model section, and validation of the VAR model). We note that unlike ordinary regression analysis, the VAR model allows for dependencies between model variables that may not be fully known 31 . For example, the number of COVID-19 cases, the degree of social distancing, and the amount of commercial activity are all interlinked as shown in Fig. 4 . If particular variables are known to be dependent (or independent), the estimation of the VAR model can be further fine tuned by constraining a sparsity structure on the regression matrix (defined in the Methods section). The procedure to identify such Figure 5 . Trace of the evolution of daily new confirmed cases and social distancing, and the associated rate of reduction in electricity consumption for two representative metropolises -New York City and Philadelphia. The bubble sizes show the reduction rate of electrical consumption (with larger bubble sizes indicating more reduction in electricity consumption). The number of COVID-19 cases and the size of the stay-at-home population are smoothed by a weekly moving average to properly extract the trends. Both cities follow a fast developing period in March 2020 and a more stable period afterwards. dependencies and impose suitable constraints on the model estimation are detailed in Supplementary Method SM-2. Further, under-reporting or inaccuracies in the reported number of COVID-19 cases are not of concern in this model, since changes in electricity consumption are only dependent on individual and policy-level responses to the observed (reported) number of COVID-19 cases (the true number of cases is unknown, and does not therefore influence changes in behaviour that contribute to a reduction in electricity consumption). We now examine the VAR model using the variance decomposition and impulse response analyses as described in Supplementary Method SM-4. Fig. 7 shows the results of the analyses for New York City, Philadelphia, and Houston (the analyses for other cities are shown in Supplementary Fig. S-5 ). Figs. 7-a,c,e show the variance decomposition of the VAR model, which indicates the contribution of different influencing factors to changes in electricity consumption. The VAR model can be further fine-tuned to enhance robustness by selecting the most significant influencing parameters, as described in Supplementary Method SM-5. We observe that decrease in commercial activity, represented by the mobility in the retail sector, is the most important factor influencing the decrease in electricity consumption since the onset of COVID-19. We also observe that the stay-at-home population can only explain a small part of the reduction in electricity consumption; one possible reason is that the change of work places from offices to homes quantified by the stay-at-home population might cause transfer of electricity consumption to the residential sector rather than a reduction. Figs. 7-b,d,f present the results of impulse response analysis on the VAR model, which describes the dynamical evolution of the reduction in electricity consumption that would result from a unit shock (1% increase or decrease) in one influencing factor. For example, in New York City, a 1% increase in the number of daily new COVID-19 cases results to a 0.29% reduction in electricity consumption in the steady state. From this analysis, the mobility in the retail sector again emerges as the most significant and robust influencing factor affecting electricity consumption across all cities (a 1% decrease in the mobility in the retail sector results in a reduction in electricity consumption ranging from 0.09% to 1.25% in the steady state.) Figs. 7-b,d,f also present important insights about the sensitivity of the electricity consumption in different cities that may not be apparent from the analysis in Table 1 . For example, Fig. 7 -f indicates that the change in electricity consumption in Houston is very sensitive to variations in the level of commercial activity, despite the magnitude of the change in electricity consumption not being very significant (Table 1) . Therefore, such cross-domain insights that are not readily available from traditional analyses may need to be considered in evaluating policy decisions pertaining to the electricity sector. Fig. 7 also indicates the "delayed" impact of various influencing factors on electricity consumption. On an average, the most pronounced impact of an increase in COVID-19 cases and social distancing is seen after a delay of one week. These findings quantify the dynamics of the interplay between the rise in the number of COVID-19 cases, increased social distancing, and reduced commercial activity, in influencing electricity consumption in the U.S. We introduced a timely open-access easy-to-use data-hub aggregating multiple data sources for tracking and analyzing the impact of COVID-19 on the U.S. electricity sector. The hub will allow researchers to conduct cross-domain analysis on the electricity sector during and after this global pandemic. We further provided the first assessment results with this data resource to quantify the intensity and dynamics of the impact of COVID-19 on the U.S. electricity sector. This research departs from conventional power system analysis by introducing new domains of data that would have a significant impact on the behavior of electricity sector in the future. Our results suggest that the U.S. electricity sector, and particularly the Northeastern region, is undergoing highly volatile changes. The change in the overall electricity consumption is also highly correlated with cross-domain factors such as the number of COVID-19 confirmed cases, the degree of social distancing, and the level of commercial activity observed in each region, suggesting that the traditional landscape of forecasting, reliability and risk assessment in the electricity sector will now need to be augmented with such cross-domain analyses in the near future. In this context, our analysis can be used as an indicator towards predicting changes in the electricity sector that will arise during the process of re-opening the economy. We also find very diverse levels of impact in different marketplaces, indicating that location-specific calibration is critically important. This work opens new directions for analysis of the electricity sector as higher quality cross-domain data sets, such as public health data on the number of COVID-19 hospitalizations or critical patients, or in-depth behavioral data on social distancing become available. In order to obtain cross-domain insights about the impact of COVID-19 on the electricity sector, we integrate data from all U.S. electricity markets with other heterogeneous data like weather, COVID-19 public health, satellite imagery, and cellular location data. The original sources for each dataset are provided in the Data and Code Availability section. Although all seven U.S. electricity markets have established websites for public information disclosure, their download centers, database structures 9/30 and user interfaces differ widely. Further, file formats, definitions, historical data availability and documentations are also extremely diverse across these markets, making it difficult to integrate this data into a unified framework. The major challenges in integrating data across different electricity marketplaces are as follows. • Some data are stored in hard-to-find pages without necessary navigation links. For example, CAISO's generation mix data are recorded in daily renewables and emissions reports, but excluded in its OASIS data platform. The download links for generation mix data in MISO can only be found by carefully checking the market report directories. • Some data are not packed and collected in an aggregated file for the requested date range. A batch downloader is needed to download these data files one by one, and then aggregate them into the desired single file. This is the case with NYISO and ISO-NE data. ISO-NE attempts to provide an interface to support data combination, but no more than 14 days of data can be pulled in one data request and the Captcha verification code system slows down the process. • Very inconsistent definitions and abbreviations are used among different markets. Sometimes, the same concepts used by different data categories don't follow the same terminology even within the same market. Some examples include: the generation fuel categories mismatch among all markets, some markets applying local time zones while others use Greenwich Mean Time, and the load zone names in ERCOT's daily files and annual archived files being different even for the same zone. • Geographical information often lacks documentation. For example, there are 17 load zones recorded in the SPP market, but their accurate locations are unknown. • The data quality is not reliable. Data redundancy, duplicate data and missing data are common problems across all markets. A typical example is LMP data. In ERCOT and ISO-NE, no price filter function is provided, and most data in large files is not needed. As shown in Fig. 1 -a, we design a processing flowchart to reorganize and harmonize all heterogeneous data sources, following three principles -data consistence, data compaction, and data cleaning, as follows. For data consistence, we design parsers to transform the data to a unified format, followed by geocoding to match the time scales of COVID-19 case load and weather data to the electricity market data, cleaning of redundant data, and file compression. Data Consistence: • For each source of electricity market files, a specific parser is designed to transform the data into a standard long table with date and hour indices. After processing by the parser, raw data from different markets is converted to a unified format. • Geocoding is adopted to match the geographical scale of three sources. • In the final labeling step, all the field names of data files are translated to the corresponding standard name from a pre-selected name list. Data Compaction: • Redundant data are dropped by parsers, and the packing step transforms the standard long tables into compact wide tables by pivoting the hour indices as new columns. Usually, the compact wide table can achieve more than 10x file compression rate compared to the un-processed raw files. • COVID-19 cases data are aggregated to the scale of market areas. • Typical ASOS stations are selected with consideration of geographical distribution and missing data rate. This step selects two to three stations for each market area, whose average missing data rates are below 0.5%. • The minute-level weather observations are re-sampled into an hourly basis to align with the resolution of market data. Data Cleaning: • Single missing data (most frequent) are filled by linear interpolation. For consecutive missing data (for example, consecutive missing dates, which are very rare), data from the EIA or EnergyOnline are carefully supplemented. • Outlier data samples are automatically detected when they are beyond 5 times or below 20% of the associated daily average value. Exceptions such as price spikes and negative prices in LMP data are carefully handled. • Duplicate data are dropped, only the first occurrence of each data sample is kept. Note that COVID-EMDA + follows the third normal form (3NF) rule to organize data in a more efficient way. This is achieved by dropping all the transitive functional dependencies in the original sources (1NF), and the data redundancy, thereby effectively reducing conflicts. The backcast model is used to estimate the counterfactual electricity consumption profile in the absence of the COVID-19 pandemic, so that the difference between a backcast model and the actual metered electricity consumption can be used to quantify the impact of the pandemic. A backcast model is expressed as a function that maps potential factors that may affect electricity consumption level, including weather variables (such as temperature, humidity and wind speed), date of year, and economic prosperity (yearly GDP growth rate) to the estimated electricity consumption. Given a group of backcast models, ensemble forecasting is widely recognized as the best approach to provide rich interval information. A group of backcast models for the daily average electricity consumption can be described bŷ whereL md is the estimated daily average electricity consumption for month m and day d,f i is the ith backcast model, T mdq , H mdq , S mdq are temperature, humidity and wind speed within the selected quantiles q, and E is the estimated GDP growth rate. We typically include 25%, 50% (average value), 75% and 100% (maximum) quantiles, and the final inputs should be decided based on the data after extensive testing. With the backcast estimations, the daily reduction in electricity consumption, r md , is calculated as follows, where T = 24 is the total number of hours in one day, and L mdt is the electricity consumption metered at time t on month m and day d. Equation (2) compares the backcast and actual electricity consumption results, and can be readily extended to interval estimations by adjusting the backcast result. We split the MISO market into three parts and the SPP market into two because these markets cover very broad areas and weather observations can vary significantly within each market area. For each market (or part of a market), we formulated a basic backcast model by considering different input settings (combination of weather variables and calendar variables with various time lag), different prepossessing techniques (normalization) and different model architectures (neural networks with various configurations, support vector machine, random forest and polynomial regression). After comparing different models, a structure with a four-layer fully-connected neural network, ReLU activation function, and L2 regularization, trained by an early stopping rule showed the best performance in terms of accuracy and robustness. Afterwards, a random search of the hidden layer numbers within ±20% fluctuation in each layer, was implemented to generate a series of backcast models. We searched 500 models for each market (or part of market) and selected the top 100 according to their estimation accuracy. Note that the mean absolute error (MAE) of these models is 2.2 − 3.4%, which is sufficiently good for one-year-ahead hourly estimation. Vector autoregression (VAR) 30 is a stochastic process model that can be used to capture the linear correlation between multiple time-series. We model the dynamics of reduction in electricity consumption using a Vector Autoregression (VAR) model of order p as follows: where . . . . in which A i is the regression matrix, x 1 stay-at-home population, median home dwell time rate, population of on-site workers, mobility in the retail sector and etc., C and E t are respectively column vectors of intercept and random errors, and the time notation t − p represents the p-th lag of the variables. Multiple pre-estimation tests are applied to check the pre-requisite properties of the input data such as stationarity and causality (see Supplementary Method SM-1). In order to avoid illogical causal relationships, some coefficients in the VAR model are restricted to be zero during the estimation (see Supplementary Methods SM-2). We traverse all combinations of input parameters and select the optimal model (see Supplementary Method SM-5). Several post-estimation tests are applied to verify the reliability of the fitted VAR model (see Supplementary Methods SM-3). Further we perform two analyses to interpret the model: (i) impulse response analysis and (ii) forecast error variance decomposition (more details in Supplementary Method SM-4). Impulse response analysis is used to describe the evolution of the VAR model's variable in reaction to a shock in one or multiple variables. Forecast error variance decomposition is to interpret the VAR model by determining the proportion of each variable's contribution to the forecast error accounting for exogenous shocks to the other variables. The COVID-EMDA + data hub and codes for all the analyses in this paper are publicly available on Github 13 . The supporting team will collect, clean, check and update the data daily, and provide necessary technical support for unexpected bugs. In the Github repository, the processed data are shared along with the original data and their corresponding parsers (written in Python). Several simple quick start examples are included to aid beginners. The original data sources for the COVID-EMDA + data hub are as follows. • Electricity market data: Data pertaining to load, generation mix and day-ahead locational marginal price (LMP) are obtained from the California (CAISO) 15 , Midcontinent (MISO) 16 , New England (ISO-NE) 17 , New York (NYISO) 18 , Pennsylvania-New Jersey-Maryland Interconnection (PJM) 19 , Southwest Power Pool (SPP) 20 , and the Electricity Reliability Council of Texas (ERCOT) 21 . Electricity market data from the Energy Information Administration (EIA) 22 and EnergyOnline company 23 is used to improve quality and fill in missing data. • Weather data: The original weather data are obtained from the Iowa State University 24 . Various kinds of weather data (temperature, relative humidity, wind speed and dew temperature) are collected from Automated Surface Observing Systems (ASOS) stations that are supported by the National Oceanic and Atmosphere Administration (NOAA), and can be extracted through an interactive website 24 . • COVID-19 public health data: The original sources of confirmed COVID-19 case numbers and deaths is the John Hopkins University dataset 25 , which contains county-level confirmed case and death numbers from January 22, 2020 (first U.S. case) onwards. • Cellular location data: The original cellular phone location dataset is derived from SafeGraph 26, 27 , a data company that aggregates anonymized GPS location data from numerous applications by census block group in order to provide insights about physical locations. The whole original dataset contains two major subdatasets: (i) social distancing metric and (ii) pattern of visits to Points of Interest (POIs). In the social distancing metric dataset, "home" is defined as the location of users at midnight, and "full-time workplace" is defined as the non-home location at which users spend more than 6 hours during daytime. The pattern dataset mainly contains (i) base information of POIs of 168 categories, including location name, address, brand association, etc., and (ii) information of daily visits and dwell time in POIs. • Satellite imagery data: In our study, the NASA VNP46A1 "Black-Marble" 28 dataset is selected as the source of satellite imagery data for its high resolution, public availability and daily update. VNP46A1 is collected by the NASA Suomi NPP sun-synchronous remote sensing satelite 33 which has an orbiting period of 101.44 minutes. This satellite measures the surface light radiation at a constant resolution of 500 meter per sample and samples daily at around local time mid-night for every location across the globe. This document presents the detailed analysis procedures and additional results to supplement the main text, and is organized into three sections as follows : 1) The following figures present visualizations of the NTL and cellular phone datasets in several hotspot cities, indicating the change of electricity consumption and mobility of the population in each city. Fig. S-8 shows the reduction in night-time light brightness, providing a visual representation of the effect of COVID-19 on electricity consumption level in major cities, as the drop in light intensity is obvious and significant. Fig. S-9 depicts a significant increase in the social distancing level indicating the change of people's mobility amidst the pandemic, with regional differences based on stringency and effectiveness of stay-at-home policies. Fig. S-11 shows the trace of the reduction in electricity consumption, new confirmed COVID-19 cases, and stay-at-home population in Boston, Houston and Kansas City, to supplement the result in Fig. 4 of the main body. Recent progress in on-board sensors and data processing algorithms for remote sensing satellites has opened up many opportunities for monitoring and analyzing human activities on the surface of the Earth and characterizing the impact of human activities on the environment, using satellite data on emission, radiation, atmosphere, vegetation, and water bodies. Among the wide range of available data, Night-Time Light (NTL) has been well recognized as a valuable and unique source of data for understanding the changes in human footprints and economic dynamics 34 . For our study, the NASA VNP46A1 "Black-Marble" 28 dataset is selected as the data source for its high resolution, public availability and daily update. VNP46A1 is collected by the NASA Suomi NPP sun-synchronous remote sensing satellite 33 which has a orbiting period of 101.44 minutes. This satellite measures the surface light radiation at a constant resolution of 500 meter per sample and samples daily at around local time mid-night for every location across the globe. This dataset has been used in power system studies from the perspectives of outage detection 35 and grid restoration 36 . The NTL dataset is used in this study as a tool for visualizing the impact of COVID-19 on electricity consumption. We conduct a comparative study of the impact of COVID-19 on artificial nightlights for representative metropolises in different RTO regions. Specifically, we focus on the cities of Boston, New-York City, Los Angles, and Houston. For each city we select a typical day in both February (before the COVID-19 outbreak) and in April (during the outbreak). The two representative snapshots selected for each city are taken from the same day-of-week and time-of-day, when the sky is clear of cloud. The raw at-sensor-Day-Night Band (DNB) data is pre-processed using the following procedures to reduce disturbances: 1. Manually locate the rectangle containing the targeted city on the tile-level NTL dataset. 2. Scan the raw data for abnormal pixels (indicated by a pixel-level Quality Flag), and approximate it by taking the average of neighboring pixels. 3. Scale every pixel using the corresponding moon illumination fraction and pixel-level lunar angles of that day to reduce disturbances from the moon. 4 . Set pixels that have extremely low light intensity (< 10 nW · cm −2 · sr −1 ) to 0, to eliminate random ambient noises. 5. Apply a 5x5 low-pass kernel filter to smooth the image. 6. Map the light intensity values of each pixel to color using a colormap and plot on axis. The processed NTL images are presented in Supplementary Fig. S-1 . The reduction in night-time light brightness provides a visual representation of the effect of COVID-19 on electricity consumption level in major cities, as the drop in light intensity is obvious and significant. The original cellular phone location dataset is obtained from SafeGraph 26, 27 , a data company that aggregates anonymized GPS location data from numerous applications by census block group in order to provide location information. The original dataset contains two major sub-datasets: (i) social distancing metrics and (ii) pattern of visits to Points of Interest (POIs). The social distancing metric dataset is generated using a panel of GPS signals from anonymous cellular phones. Note that "home" is defined as the common nighttime location of each cellular phone over a 6 week period, and part-time and full-time workplaces are defined as the non-home locations where users spend from 3 − 6 hours and ≥ 6 hours respectively between 8 am and 6 pm in local time. The following are the features selected for our analysis: • Basic Information: (i) unique 12-digit FIPS code for the Census Block Group; (ii) start and end time for the measurement period (namely 24 hours); (iii) count of devices whose homes are in the Census Block Group. • Completely Stay-at-home Device Count: the number of devices that never leave "home" during the measurement period, out of the total count of devices in the Census Block Group. • Median Home Dwell Time: the median dwell time at "home" in minutes during the measurement period, for all devices in the Census Block Group. • Part-time Work-on-Site Device Count: the number of devices that go to part-time workplaces during the measurement period, out of the total count of devices in the Census Block Group. • Full-time Work-on-Site Device Count: the number of devices that go to full-time workplaces during the measurement period, out of the total count of devices in the Census Block Group. We aggregate the daily social distancing data by county. Denote the total count of the completely stay-at-home devices in a county as C 1 , the median value of the median home dwell time in a county as C 2 , the total count of the part-time work-on-site devices in a county as C 3 , the total count of the full-time work-on-site devices in a county as C 4 and the total count of devices in a county as C. Then, we define C 1 /C, C 2 /1440, C 3 /C, and C 4 /C as the county-level "completely stay-at-home rate", "home dwell time rate", "part-time work-on-site rate" and "full-time work-on-site rate" respectively, which will be used for the VAR model. The pattern dataset is a place traffic and demographic data aggregation available for about 4 million POIs that contains the frequency of visits to various POIs, the dwell time, the residence location of visitors, etc. The following are the features selected for our analysis: • Basic Information: (i) the unique and consistent ID tied to each POI; (ii) name of the POI; (iii) physical address; (iv) postal code; (v) brand; (vi) start and end time for measurement period (about one week) in local time; (vii) associated category of the POI. • Visits by Day: the number of visits to the POI each day over the covered time period. The daily pattern data for "retail" POIs are collected as defined in Supplementary Note SN-3. We aggregate and sum up the total number of daily visits by county, which will be used as "retail mobility" data in the VAR model. After verifying the input variable time-series using the statistical tests described in Supplementary Method SM-1, we have a set of de-trended stationary input time-series that do not have long-term correlation among them. We model the dynamics of load reduction as a Vector Autoregression (VAR) model of order p as follows: where . . . . in which C and E t are respectively column vectors of intercept and random errors, x 1 t represents the target output variable at time t, namely the load reduction amount we wish to model, x 2 t , ..., x n t represent the selected n − 1 parameter variables including the number of COVID-19 cases, completely stay-at-home rate, median home dwell time rate, etc, and the time notation t − p represents the p-th lag of the variables. The coefficients a k i, j are calculated separately for each variable using the Ordinary Least Square (OLS) estimator: Then, we can aggregate and concatenate all a k i, j to obtain the regression matrix A k , 0 ≤ k ≤ p. If the result from the Granger causality test suggests that there may exist undesirable causal relationships between variables, we can impose constraints on the OLS to eliminate such relationships from appearing in the VAR model. For example, if the Granger test suggests that variations in load reduction could be a causation of change in the stay-at-home population, which is intuitively illogical, we can restrict the corresponding entries of the VAR model coefficient matrices that describe such a relationship to zero during the OLS computation. With this Restricted VAR model, we ensure that the final model does not include unwanted relationships between parameter variables. To verify the restricted VAR model, we need to guarantee stationarity, non-autocorrelation, and normality of the residual data, which is defined as where e t and x t are the residuals and observation vectors respectively on day t, k is the user-defined maximum lag, namely the order of the VAR model, and {Â i } k i=1 are the estimated coefficient matrices derived from the estimated VAR model. In addition to verifying the stationarity of the input time-series data in Supplementary Method SM-1, the ADF test is also used to check if the residual data are non-stationary and possess a unit root. The endogeneity of the residual also needs to be verified, since the existence of endogeneity may render the regression result untrustworthy. Therefore, the Ljung-Box test 39, 40 is used to test whether any of a group of autocorrelations of the residual time-series are different from zero. In this test, the null hypothesis is that the data are independently distributed while the alternative hypothesis is that the data are not independently distributed and exhibit serial correlation. The test statistic 40 is where n is the number of residual data, h is the number of lags selected to be 40, and ρ i is the sample autocorrelation value at lag i. With significance level 5%, the critical region for rejection of the null hypothesis is Q > χ 2 0.95,h . The Durbin-Watson statistic ? is another test statistic used to detect the presence of autocorrelation at lag 1 in the residuals of the VAR model. The null hypothesis is that the residuals are serially uncorrelated while the alternative is that they come from a first order autoregression process. Normality tests are used to determine if a dataset is well-modeled by a normal distribution. The null hypothesis is that the data comes from a normal distribution. Specifically, the Statsmodels module implements the normality test of the VAR residuals by the Jarque-Bera (JB) test ? . The JB statistic has a χ 2 distribution with two degrees of freedom if the data comes from a normal distribution. Therefore, the null hypothesis of the JB test is the joint hypothesis that the skewness is 0 and the excess kurtosis is 0. With a VAR model that passes all statistical tests, we need to analyze and interpret the results of the model. We perform the following analyses on the VAR model. Impulse response analysis ? is an important step to describe the evolution of a VAR model's variable in reaction to a shock in one or more variables. For the p-th order restricted VAR model of the form where x t , c and e t are n dimensional column vectors, A i is p × p dimensional matrix. In the impulse response analysis, we set the impulse response function R(t) as with R(0) = [0, . . . , 1, . . . , 0] T in which only one user-defined element equals to 1. It is particularly useful to consider the impulse response function as a predictive indicator that can forecast the dynamic behavior of electricity consumption in response to a change in any exogenous influencing factor, given a certain initial state. For example, suppose we have another model describing how public health policies impact social mobility, we can further combine these two models and simulate the impacts of public policies on the electricity consumption. Forecast error variance decomposition (FEVD) ? is used to aid in the interpretation of the estimated restricted VAR model by determining the proportion of each variable's forecast error variance that is contributed by shocks to the other variables. For the p-th order restricted VAR model in the form where x t , c and e t are n dimensional column vectors, A i is p × p dimensional matrix. Then, it can be reformulated as where x t x t−1 x t−2 . . . and I n is a n × n identity matrix, A is a np × np matrix, and X t , C and E t are np dimensional column vectors. The mean square error of the h-step forecast of the i-th variable can be calculated as where Σ e is the covariance matrix of e t , Φ j = JA j J T , J = [I n , 0, . . . , 0] k×kp , A j is the j-th order power of A matrix, and (i, i) denotes the i-th diagonal element of the matrix. Further, for the forecast error variance decomposition, w i j (h) is defined to represent the proportion of forecast error variance of the i-th variable at the h-th step accounted for by the shock to the j-th variable, that is, where e i is the i-th columns of the identity matrix I n×n , B k = Φ k P, and P is a lower triangular matrix in the Cholesky factorization of Σ e such that Σ e = PP T . The aforementioned methods compute the coefficients of a VAR model that is statistically robust. However, there is no explicit rule for selecting the parameter variables and the range of the training data. A numerical search on the parameter space is performed for each city to determine the optimal parameters. A list of parameters we use for iterative search is listed below. • Time-series from EMDA dataset used as input variables: [x 1 , ..., x n ] • Date range of training data: [T 1 , T 2 ] • Order of the VAR model: p • Rule to determine whether to set a restricted VAR coefficient to 0: r For each combination of possible parameters, the procedures described in Supplementary Methods SM-1,2,3,4 are performed to examine feasibility and quantify the numerical performance. After searching the parameter space, we determine the optimal combination of parameters according to the FEVD and impulse response analysis. An ideal model should be able to to explain a large proportion of the variance in load reduction, while at the same time ensuring that the signs of the trends (increase/decrease) of all impulse responses are as desired. The complete process of searching for the optimal model is presented in Algorithm 1: Utilities beginning to see the load impacts of COVID-19 as economic shutdown widens Electricity demand in the time of COVID-19 COVID-19: Potential impacts on the electric power sector Tracking the impact of coronavirus on the us power sector Impacts and implications of covid-19 for the energy industry Electric Power Research Institute. Demand impacts and operational and control center practices Sharing knowledge on electrical energy industry's first response to covid-19 Impact of covid-19 on power system operation planning COVID-19 load impact analysis (ERCOT Midcontinent Independent System Operator Reliability Subcommittee. COVID-19 impact to load and outage coordination (MISO Impacts of COVID-19 on NYISO demand Electric power markets: national overview California ISO open access same-time information system site Midcontinent Independent System Operator Market and operations: ISO express Electric Reliability Council of Texas Energy Information Administration. Hourly electric grid monitor Iowa environmental mesonet: ASOS-AWOS-METAR data download An interactive web-based dashboard to track COVID-19 in real time Black marble user guide version 1.0 (NASA Carto dark "dark matter Macroeconomics and reality Vector autoregressions A vector autoregression weather model for electricity supply and demand modeling National Aeronautics and Space Administration Applications of satellite remote sensing of nighttime light observations: Advances, challenges, and perspectives Monitoring disaster-related power outages using nasa black marble nighttime light product (ISPRS Satellite-based assessment of electricity restoration efforts in puerto rico after hurricane maria Distribution of the estimators for autoregressive time series with a unit root Information theory and an extension of the maximum likelihood principle Distribution of residual autocorrelations in autoregressive-integrated moving average time series models On a measure of lack of fit in time series models Automotive Parts, Accessories, and Tire Stores Goods Stores 17. Lawn and Garden Equipment and Supplies Stores 18 Sporting Goods, Hobby, and Musical Instrument Stores Utilities beginning to see the load impacts of COVID-19 as economic shutdown widens Electricity demand in the time of COVID-19 COVID-19: Potential impacts on the electric power sector Tracking the impact of coronavirus on the us power sector Impacts and implications of covid-19 for the energy industry Electric Power Research Institute. Demand impacts and operational and control center practices Sharing knowledge on electrical energy industry's first response to covid-19 Impact of covid-19 on power system operation planning COVID-19 load impact analysis (ERCOT Midcontinent Independent System Operator Reliability Subcommittee. COVID-19 impact to load and outage coordination (MISO Impacts of COVID-19 on NYISO demand Electric power markets: national overview California ISO open access same-time information system site Midcontinent Independent System Operator Market and operations: ISO express Electric Reliability Council of Texas Energy Information Administration. Hourly electric grid monitor Iowa environmental mesonet: ASOS-AWOS-METAR data download An interactive web-based dashboard to track COVID-19 in real time Black marble user guide version 1.0 (NASA Carto dark "dark matter Macroeconomics and reality Vector autoregressions A vector autoregression weather model for electricity supply and demand modeling National Aeronautics and Space Administration Applications of satellite remote sensing of nighttime light observations: Advances, challenges, and perspectives Monitoring disaster-related power outages using nasa black marble nighttime light product (ISPRS Satellite-based assessment of electricity restoration efforts in puerto rico after hurricane maria Distribution of the estimators for autoregressive time series with a unit root Information theory and an extension of the maximum likelihood principle Distribution of residual autocorrelations in autoregressive-integrated moving average time series models On a measure of lack of fit in time series models The "retail" data refers to aggregation of several categories of POIs from the "pattern of visits to POIs" dataset 27 . We select the following 25 categories among a total of 168 POIs as indicators of mobility in the retail sector: The Statsmodels ? module in Python is used to implement the estimation of coefficients in the restricted Vector Autoregression (VAR) model and the corresponding statistical tests and analyses. We describe the key steps in estimating and verifying the VAR model as follows. Collecting electricity market data, weather data, COVID-19 public health data, satellite image data and cellular location data, we pre-process the following variables as input candidates for the restricted VAR model:• Load Reduction: Logarithm of the amount of load reduction (MW) -samples with negative reduction (increase) are dropped for consistency.• New Daily Confirmed Cases: Logarithm of the original count.• Stay-at-home Population: Logarithm of the number of devices that stay at home completely. • Mobility in the Retail Sector: Logarithm of the number of visitors to retail locations as defined in Supplementary Note SN-3. The stationarity of time-series data is a prerequisite to calibrate a vector autoregression (VAR) model. Therefore, the Augmented Dickey-Fuller (ADF) test 37 , a commonly used unit root test, is used to test whether a time-series variable is non-stationary and possesses a unit root. The Akaike Information Criterion (AIC) 38 is selected to test the null hypothesis that a unit root is present in a time-series sample. The Augmented Dickey-Fuller test is carried out for each of the multiple time-series data that are candidate variables for VAR model estimation. Each time-series is differenced to improve stationarity. The cointegration test ? is used to test the long-term correlation between multiple non-stationary time-series. In addition to the ADF test which tests the stationarity of each de-trended input time-series, a cointegration test is used to detect potential long-term correlation among the original inputs, which is signified by the presence of cointergation. If cointegration is detected, the selected tuple of input time-series is not appropriate for VAR modelling and needs to be dropped. Granger causality ? is a probabilistic method to estimate the causality relationship among two variables represented as time-series.The key intuition behind the Granger Causality test is the assumption that future events cannot have causal effects on the past. To study the causality relationship between two variables, x and y, we test if the lagged series x t−n , n ∈ Z + will affect the current value y t for each time-step t, as concurrent and future values x t+n , n ∈ Z cannot affect y t . The Granger Causality test is used to estimate the causal relationship between selected time-series variables shifted by the appropriate lag values. If significant counter-logical casual relationships are identified, we impose constraints on VAR modeling such that the causalities are eliminated in the final VAR model.