1

A Sparse Linear Model and Significance Test
for Individual Consumption Prediction
arXiv:1511.01853v3 [stat.ML] 21 Feb 2017

Pan Li, Baosen Zhang, Yang Weng, and Ram Rajagopal

Abstract
Accurate prediction of user consumption is a key part not only in understanding consumer flexibility
and behavior patterns, but in the design of robust and efficient energy saving programs as well.
Existing prediction methods usually have high relative errors that can be larger than 30% and have
difficulties accounting for heterogeneity between individual users. In this paper, we propose a method
to improve prediction accuracy of individual users by adaptively exploring sparsity in historical data and
leveraging predictive relationship between different users. Sparsity is captured by popular least absolute
shrinkage and selection estimator, while user selection is formulated as an optimal hypothesis testing
problem and solved via a covariance test. Using real world data from PG&E, we provide extensive
simulation validation of the proposed method against well-known techniques such as support vector
machine, principle component analysis combined with linear regression, and random forest. The results
demonstrate that our proposed methods are operationally efficient because of linear nature, and achieve
optimal prediction performance.

Index Terms
Load forecasting, least absolute shrinkage and selection, sparse autoregressive model, significance
test

N OMENCLATURE
yt

Consumption at time t.

Pan Li and Baosen Zhang are with the Department of Electrical Engineering, University of Washington, Seattle, WA, 98195,
(e-mail: {pli69, zhangbao}@uw.edu).
Yang Weng and Ram Rajagopal are with the Civil and Environmental Department, Stanford University, Stanford, CA, 94035,
(e-mail: {yangweng, ramr}@stanford.edu).

2

ŷt

Estimated consumption at time t.

βi

Parameters in the regression model.

β

Parameters as a vector.

εt

Noise at time t.

y

Consumption data as a vector.

X

Regressor matrix in the regression model.

β OLS
β̂

OLS estimate of β .

β LASSO LASSO estimate of β .
β̂
λ

Tuning parameter in LASSO.

e

Fitted Residual after LASSO.

ξ

Regressor matrix in the regression model for residual.

α

Parameters in the regression model for residual.

ψ

Noise at time t in the regression model for residual.

α (λk )
α̂

Estimated α at tuning parameter λ = λk .

σ2

Variance of the noise.

σ̂ 2

Estimated variance of the noise.

F1

Test statistic.

T1

Test statistic.

Fdi ,d j

F distribution with di and d j as parameter.

A

Active set.
I. I NTRODUCTION

Electric load forecasting is an important problem in the power engineering industry and
have received extensive attention from both industry and academia over the last century. Many
different forecasting techniques have been developed during this time. The authors in [1] present a
comprehensive literature review on different methods related to load forecasting, from regression
models to expert systems. Time series methods are further discussed in [2]. A thorough research
on load and price forecasting is presented in [3]. A common theme among many of the established
methods is that they are used to forecast relative large loads, from substations serving megawatts
to transmission networks serving more than gigawatts of power [4]. As the grid shifts to a more
distributed system, the need of accurate forecasting for smaller sized loads is becoming increasing

3

important.
Recent advances in technology such as smart meters, bi-directional communication capabilities
and distributed energy resources have made individual households active participants in the power
system. Many applications and programs based on these new technologies require estimating the
future load of individual homes. For example, state estimation algorithms for distribution systems
require pseudo-measurements [5], and these are provided by load forecasts. Another important
class of application is demand response and dynamic pricing programs, where users’ future
demand are needed to design appropriate incentives [6], [7], [8], [9], [10].
In contrast to large aggregated loads, individual load forecasting is less developed. The current
state of forecasting algorithms falls under three broad classes: simple averaging, statistical
regression methods, and Artificial Intelligence (AI). They are listed in increasing order of
prediction accuracy and decreasing in order of model simplicity. Simple averaging is intuitively
pleasing since it is based on the mean of the previous similar days, but is often not very accurate.
On the other end of the spectrum, AI methods can be extremely accurate, but it is difficult to
associate the obtained parameter with the input data, i.e., past consumption. In this paper, we
propose an algorithm that can achieve the performance of the state of the art AI methods, but
retains the simplicity of linear nature of regression methods.
The algorithm we propose is based on the well known Least Absolute Shrinkage and Selection
(LASSO) algorithm in statistics and signal processing [11]. It shrinks parameter estimates towards
zero in order to avoid overfitting as well as pick up the most relevant regressors. Operationally,
it minimizes the usual sum of squared errors, with a bound on the sum of the absolute values of
the coefficients. LASSO is preferable in a setting where the dimension of the features are much
higher than the size of the training set [12], [13]. Furthermore, it can be easily extended to suit
for more scenarios, i.e., group LASSO or adaptive LASSO [14], [15]. Due to its nice statistical
properties and its efficiency, LASSO is applied to many different disciplines and has shown its
superiority in sparsity recovery [16], [17].
In particular, we formulate the load forecasting problem as learning parameters of a sparse
autoregressive model using LASSO. This sparse autoregressive model automatically selects the
best recurrent pattern in historical data by shrinking irrelevant coefficients to zero. By selecting
the correct features, the algorithm improves the order selection in autoregressive models, and
as we will show using real load data from Pacific Gas and Electric Company (PG&E), it also

4

greatly improves the prediction accuracy of current regression models. For example, taking the
Mean Absolute Percentage Error (MAPE) as a metric, autoregressive model with lag order one
(AR(1)) has a MAPE of 33.9 %, based on a pool of 150 users. The proposed method reduces the
MAPE to 22.5 %. Nonparametric methods such as using the consumption data at the same time
during last week (LW method) are more intuitive but their performance are not robust subject
to noise. They yield a relative error nearly up to 35% based on a pool of 50 users in our test
dataset. In addition, based on our dataset, linear regression model with the reduced feature space
from Principle Component Analysis (PCA) performs rather poorly, with a relative error greater
than 100 %, which is even worse than simply predicting zero. Overall, the proposed model with
LASSO achieves a much smaller error and maintains the intuitive nature of linear methods.
In addition, we observe that given the set of coefficients, utilities and planners can easily
identify and attribute the impact of different features on load consumption based on LASSO.
For example, suppose for a given user, our method concludes that its load at time t is mostly
determined by the past load at times t − 1, t − 3 and t − 24. This shows that this user has some
short term behavior in the one hour range, medium term behavior in the 3 hour range, and a
daily cycle partner that repeats every 24 hours. This information can be interpreted by utilities to
better understand this user’s consumption patterns and potentially identify appropriate demand
programs.
Apart from LASSO, data from other users can be leveraged for forecasting in the proposed
framework. Intuitively, this means that knowing the past history of user j improves forecasting
of the user i. For a given user, we use a sequential hypothesis testing procedure to find the
best other user’s historical data to include in the algorithm. We give a rigorous derivation of the
hypothesis testing procedure and quantify the confidence of including other users. This allows
us to show that the procedure is optimal, in the sense that given user i, it finds user j, whose
historical data improves load forecasting for user i the most among all users.
We derive rigorous theoretical justification for our methods as well as provide extensive
simulation studies with respect to several well studied prediction methods. In particular, we
compare against AutoRegression (AR), Exponential Smoothing (ES), Support Vector Machine
(SVM), linear regression with PCA, Random Forest (RF) and Neural Network (NN) model. We
also include two more nonparametric methods, i.e., using ten previous days’ average consumption
as prediction and using last week’s (LW as mentioned earlier) consumption data as prediction.

5

Using a user’s own historical data, our proposed method and RF both reduce prediction error
by 30% compared to other predictors. Our proposed method is simpler in nature than RF since
the latter is a generic machine learning technique that relies on a random ensemble of decision
trees [18], [19]. Therefore, our method is useful to system operators in policy decisions without
sacrificing prediction accuracy. By adding the historical data of another user, we can further
improve the prediction accuracy.
A. Contribution
The contribution of this paper is two fold.
•

First, we apply LASSO in consumption prediction. LASSO has been studied in literature to
jointly predicts price, load and renewables together [20]. However, there is little discussion
on how LASSO can be applied to consumption forecasting without using other side information. Our paper has a thorough discussion on how LASSO is used in obtaining a sparse
autoregressive model and compares its performances with several other popular prediction
methods in literature. The simulation results show that LASSO achieves competitive prediction performance compared to nonlinear machine learning algorithms such as RF or NN.
In addition, LASSO is computationally much faster and easier to be understood by human
operators.

•

Second, we propose a significance test that can leverage other users’ consumption data for
prediction. This testing procedure differs significantly from standard clustering algorithms
since it looks for the most “predictive” user, not necessarily the most similar user.

The rest of the paper is organized as follows. Section I-B analyzes related work in short term
load forecasting. Section II presents the autoregressive model for time series analysis. Section III
introduces LASSO type linear regression model. Section IV proceeds with the significance test
to pair users in order to improve prediction performance. It describes the significance test for
LASSO, i.e., covariance test, to select the most significant user to form the pair with the current
user. Section V introduces the evaluation criterion for prediction and details the simulation of the
proposed methods compared to several other popular prediction methods. And finally Section
VI concludes the paper and draws avenues for future work.

6

B. Literature review
There exists an extensive literature on short term load forecasting [21], [1]. In summary,
research on energy consumption prediction can be divided into three groups [1], including simple
averaging models, statistical models and AI models.
The simplest approach is to employ moving average [22]. Such models make predictions
on mean of consumption data from previous similar days [23]. AI type methods (e.g. NN or
RF) yield high accuracy at the cost of complexity of the system, which may lead to overfitting
[24]. Other drawbacks include difficult parametrization and non-obvious selection of variables.
Statistical methods sit in between the previous methods in terms of complexity and accuracy,
and include regression models, probabilistic approach applied to regression models, and time
series analysis such as autoregressive models.
In statistical methods, regression models combine several features to form normally a linear
function. In [25], the authors build a regression tree model with weather data to predict consumption. SVM is used in [26]. Gaussian process framework for prediction mitigating the uncertainty
problem is proposed recently in [27].
Advanced AI models have also been used to facilitate load forecasting. Authors in [28] propose
a NN model to handle uncertainties associated with forecasts. Similarly, authors in [29] study
the interruption prediction with weather conditions by simulating a NN. Besides these work,
the authors in [30] explore kernel-based learning techniques to forecast both commercial and
industrial load. In addition, feature selection techniques are applied in [31] to predict electricity
load and prices.
Apart from these advanced machine learning tools, another type model involves time series
based methods [32]. An overview can be found in [33]. Models such as ES [34] and the
Autoregressive Integrated Moving Average (ARIMA) model [35], [36] are all time series based.
ARIMA is a widely used time series based prediction methods and have been adopted in [37],
[38]. In these papers the authors apply ARIMA to predict either the electricity price or the
electric vehicle demand. Apart from these work, authors in [20] propose a vector autoregressive
model to include renewables, prices and loads together with sparsity recovery. In this work, the
authors have also explored LASSO to obtain a sparse linear model. In addition, to extend from
linearity to nonlinearity, [39] addresses a mixed model combining ARIMA model to deal with

7

the linear part and neural network with the nonlinear one. In our work, we recover the sparsity
for univariate time series and multivariate time series under the framework of autoregressive
models.
II. AUTOREGRESSIVE

MODEL

Autoregressive models are widely used for prediction and inference of time series data. Here
we adopt a linear autoregressive model of the hourly consumption of a single user, where future
demands yt depend linearly on historical data yt−i plus random noise:
I

yt = β0 + ∑ βi yt−i + εt .

(1)

i=1

In this model, yt denotes the demand of the user at time t, βi is the coefficient for order (lag
order) i in the autoregressive model (it represents the weight of each historical demand data yt−1
in predicting future demand yt ), and εt is an additive random Gaussian noise. The time index t is
measured in hours and the noise is identically and independently distributed at different hours.
Note that in this paper we denote time series data by notation {•t }, where subscript t refers to
the time slot in this time series data. In addition, I is the number of orders that we include in
the model. An autoregressive model with maximum order I is denoted by AR(I).
To use the model in (1) for prediction, the standard approach is to use Ordinary Least Squares
(OLS) to estimate the coefficients βi , i ∈ I. By convention, we write [yt−1

yt−2

···

yt−I ]T

as a vector denoted by xt . Using this notation, the model in (1) is written in a compact matrix
form:
y = Xβ + ε ,
where y = [yt
= [εt

εt+1

yt+1

· · · ]T , X is a matrix where t th row is [1 xTt ], β = [β0

(2)

β1

· · · ]T , and ε

· · · ]T . Vectors y, β and ε have dimension T . Matrix X has P columns, which we

refer to as the dimension of X.
Applying standard OLS to (2), the estimate of β is given by:

β OLS = arg min k(y − Xβ )k22 .
β̂

(3)

β

β OLS is a consistent estimator for the
Under some standard assumptions, the OLS estimator β̂
true β , meaning that the expected difference between the estimator and the true value approaches
zero when sample size becomes large [40]. This means that the bias goes to zero as the sample

8

size becomes big. The other asymptotic analysis on the OLS also applies in this case, such
as asymptotic distribution for gaussianality of estimators and significance tests; however, in
autoregressive models the estimators are typically not unbiased.

β OLS from OLS, we separate the dataset into a training
To learn and evaluate the estimator β̂
β OLS is learned from the data from the training set and the
set and a test set. The estimator β̂
estimation error is evaluated on the test set. Note that estimators may exhibit extremely good
fits on the training set but poor estimation performance on the test set, as according to [41].
III. S PARSITY

IN

AUTOREGRESSIVE M ODELS

Since the objective of OLS estimators is to minimize the sum of squared errors in the training
set, OLS achieves optimal in-sample performance. This means that adding more regressors into
(2) can always decrease the sum of squared error and better fit the data within the training set.
However, when we include too many irrelevant regressors, i.e. when we include too many lag
orders from the historical data in (2) , we are misled by the reduced in-sample bias. We will
then ignore the high variance introduced by the estimator which leads to model overfitting.
We use the PG&E dataset as an illustrative example. In this particular dataset, hourly consumption data for single households is recorded. If we use an AR(5) model for a particular household,
it will result in an average in-sample squared error of 0.0159, with an average out-of-sample
error of 0.0216, whereas AR(1) model has an average in-sample squared error of 0.0172, together
with an average out-of-sample error of 0.0174. Thus AR(1) gives better out-of-sample fitting
results. If the potential lag orders are up to 10 days, i.e., 240, then an AR(240) model would
produce large out-of-sample errors. Overall, we need to select the lag orders carefully to avoid
model overfitting. The tradeoff between out-of-sample error and the order of the autoregressive
model is shown in Fig. 1.
Determining the correct lag orders is not trivial because it is a combinatorial problem, which
is NP-hard. To this purpose, we use LASSO, which is a convex relaxation of such combinatorial
problems, to select relevant lag orders. The intuition for using LASSO is to get a sparse
autoregressive models with high orders.
More formally, LASSO is a shrinkage and selection method for linear regression described
in (2) [11]. It shrinks parameter estimates towards zero in order to avoid overfitting as well
as picks up the most relevant regressors. Operationally, it minimizes the usual sum of squared

9

0.026
In-sample
Out-of-sample

0.024

Error

0.022
0.02
0.018
0.016
0.014

1

2

3

4

5

Order in autoregressive model
Fig. 1: Error comparison between different autoregressive models.

errors, with a bound on the sum of the absolute values of the coefficients:
1
β LASSO = arg min k(y − xβ )k22 + λ kβ k1 ,
β̂
2
β

(4)

where λ is a tuning parameter to control the level of sparsity in the solution. In essence, it
controls the number of past consumptions included in the prediction of future consumptions

β LASSO is, which means that less
(i.e. sparsity). The bigger λ is, the more sparse the solution β̂
number of data points from the past is used to predict the future. This parameter balances the
training and testing performances of the model. When λ = 0, the solution is the same as in (3)
as a traditional autoregressive model. For practical purposes, we are using k-fold cross validation
(CV) to determine the value of λ in our simulations, where k is either 5 or 10 [41].
LASSO has gained wide spread popularity in signal processing and statistical learning, see
[42], [43], [44]. LASSO has also been applied to forecast electricity price [20], [45], but
its application to load forecasting is still a new topic. In [46], LASSO has been applied to
forecast short and middle term load. In our paper, we adopt LASSO to predict hour-ahead
consumption and provide a comparison between LASSO and some other well-used prediction
methods in literature. As discussed later in Section V, LASSO is the most preferable model when
considering model simplicity and prediction performances. Despite the fact that RF achieves the

10

best prediction performance with an average relative error less than 20%, it is a highly non
linear model and requires to tune many hyper-parameters. This is computationally expensive
and not linear in nature. On the other hand, LASSO achieves similar performance as RF and
outperforms all other linear models considered in this paper, i.e., SVM with a linear kernel,
linear regression with PCA, AR(1) and simple averaging method. It reduces the relative error
by 30% compared to these models. Thus LASSO is considered the best in terms of simplicity
and prediction accuracy based on our dataset.
Another advantage of using LASSO to recover sparsity is that it has some nice properties as

β has the same support as β and the sign of each
to sign consistency [12], [13]. This means that β̂
element in this support is recovered correctly. Therefore LASSO recovers the exact sparsity of
the underlying model. In our simulation, LASSO selects both the most recent lag orders and lag
orders with intervals of roughly 24 hours, which performs as a combination of simple averaging
and AR(·). Furthermore, LASSO also gives more intuitive results with respect to selected orders.
In our simulation, for one electricity user as an example, LASSO selects lag orders as 1, 2, 5,
6, 16, 22, 23, 24, 48, 143, 144, 160, 191, 216, 238, 240. From these orders we can observe a
clear behavior pattern of an interval of 24 hours. Some are multiples of 24 and some are not
but close to multiples of 24. We thus observe that not every lag order that LASSO picks is a
multiple of 24, otherwise we would directly employ simple averaging rather than LASSO, so
LASSO is more adaptive and flexible than simple averaging or AR(·). This implies that user
behavior at current hour depends on similar hours happened in previous days. Unlike simple
averaging which fix the lag orders at 24, 48, 72, etc., LASSO will automatically select these
orders for each individual based on their respective historical data, instead of directly imposing
fixed orders.
As can be seen from Section V, LASSO applied to autoregressive model achieves the best
prediction performance among all the linear models considered in this paper, with a relative
error as small as 22.5 %. It improves the prediction performance by 33.6 % as compared to an
AR(1) model.
IV. U SER PAIRING

BY

S IGNIFICANCE T EST

So far we have considered using historical data of an individual user for its own prediction. One
way to leverage the fact that we have many users’ data is to improve the univariate autoregressive

11

model by including other user’s historical data into the model. One popular way to perform this
is to conduct Vector Autoregression (VAR), which extends the univariate autoregressive model
to joint prediction for a vector of time series data.
To perform a complete VAR, we need to include all potentially relevant users into the
autoregressive model, which will reduce the bias but increase the variance for estimators. This
causes the same overfitting problem as occurred in univariate autoregressive model, when AR(3)
yields a worse prediction on the test set but a better fit on the training set compared to AR(1).
One possible way to overcome this problem is to first cluster similar users together and then
perform VAR for each cluster. However, consider a scenario where two time series have the
exact same values for each time slot. Then these two time series are clustered together since they
are identical. In this case, clustering them together and performance multivariate autoregressive
model does not help to enhance prediction because knowing the history of one time series would
not help predict the future values of another time series. This problem distinguishes similarity
from prediction performance, which is the focus of this section.
In this paper we focus on selecting the most relevant user to enhance prediction performance
after doing univariate LASSO selection. To this end, we adopt LASSO significance test to select
this most relevant user. In LASSO significance test, the inclusion of a particular user is based on
how well this particular user’s history data explains the fitted residual after performing LASSO
to one user’s univariate autoregressive model. This is a hypothesis test of an exponential random
variable [47]. We will discuss more details of a LASSO solution to a regression model and
the implementation of the LASSO significance test in the rest of this section. Overall from the
simulation results presented in section V, using significance test on top of univariate LASSOtype regression model improves the relative prediction error from 22.5% to 20.9%, which is
almost as good as the prediction results from RF.
A. Linear Regression Model for the Fitted Residual and LASSO Path
In LASSO significance test, we want to test if the fitted residual from LASSO solution in (4) is
indeed noise or if it can be better explained by other user’s historical data. Intuitively, We need
to test inclusion of each user’s historical data. We therefore generate LASSO autoregressive
coefficients for univariate time series in (1) and compute the residual {et } (true value minus
fitted value) for each individual user. Then we use et as the response variable in a new linear

12

regression model. To fit this new linear regression model, we include all users’ historical data at
lag order one except for this user as the regressors. In this way we will have a high dimensional
regressor matrix. Its dimension is the number of users we want to potentially include, which is
the total number of users minus one. Since we only want to include the most significant user as
the regressor, LASSO is therefore performed and furthermore tested on this new linear regression
problem. In the following, we first illustrate more details on formulating the LASSO regression
model to the fitted residual from univariate model. Then we provide more explanations on the
LASSO path with a varying λ . In the next subsection we introduce the covariance statistic and
its asymptotic distribution for hypothesis tests along the LASSO path.
Mathematically, the new linear regression problem for the residual is formulated as follows:
e = ξ α +ψ,
where e = [et,s

et+1,s

T β̂
β LASSO
· · · ]T = [yt,s − Xt,s

for the current user s, Xt,s = [yt−1,s

yt−2,s

···

(5)

T
β LASSO
β̂
yt+1,s − Xt+1,s

· · · ]T is the residual

yt−I,s ]T and yt,s denotes the consumption data

for user s at time t. The parameter α is a vector of weights that denotes how helpful each
other user’s last demand data is towards the prediction of future demand of user s. Vector
e has a length of N. Moreover, matrix ξ is the regressor matrix made by lag order one
historical data from all the other users except for the current user s. The t th row of ξ is
[yt−1,1

yt−1,2 · · · yt−1,s−1

yt−1,s+1 · · · ]. So ξ has a column length of P, which is the number

of users to be included, and a row length of N, which is the same as that of e. One user is
represented by one column of ξ and is regarded as one regressor variable. In addition, ψ is the
white noise vector with variance σ 2 . The parameter α = [α0

α1

···

αP ]T is the decision

variable for the regression problem.
We again apply LASSO for (5), to avoid overfitting by including too many irrelevant users.
Building on the discussion in Section III, we define the LASSO path as the revolution of the

α in terms of a sequence of λk ’s. The LASSO path α̂ (λk ) is given by:
estimator α̂
1
α (λk ) = arg min k(e − ξ α )k22 + λk kα k1 ,
α̂
2
α

(6)

where λk is called the knot along the LASSO path.
For different values of λk , we obtain different solutions and sparsity at different levels. The
active set at one particular value of λk is the set of all non zero coefficients estimated at that

13

α (λk ) is continuous and piecewise
value, i.e., A = {α̂ p 6= 0, λ = λk , p = 0, 1, · · · , P}. The path α̂
linear with knots at these values λ1 > λ2 > · · · > 0 [48]. With the formulation in (5), the goal
is to test if an inclusion of one user’s historical data is helpful for prediction. Mathematically
speaking, we want to test if the variables that sequentially enter the active set are statistically
significant.
B. Covariance Test
Significance test applied to the LASSO path is fundamentally different from the standard
hypothesis tests in linear regression. In standard testings with linear regression model, we
fix a regressor variable and test the significance of its inclusion into the model. This can be
done by t-test. However, the scenario is not the same in the case where we want to test the
significance of variables along the LASSO path since the variables entering the active set are
not known in advance. This means that we do not know which regressors to fix in order to
conduct t-test. This problem is addressed by the authors in [47]. Instead of standard t-test, they
propose a significance test for regressors progressively selected by the LASSO path, using the so
called covariance statistic. They have shown that the covariance statistic asymptotically follows
exponential distribution and can be used to test the significance of the entry of variables into
the active set.
Using the covariance test defined in [47], we test whether it is significant to include one
variable (one other user) into the active set against including no variables into the active set.
Specifically, we want to analyze whether the residual, i.e., e is white noise, or it can be explained
significantly by the lag one order historical data from other users.
The null hypothesis is that the true active set is an empty set, i.e., no other users can help
better predict the current user:
H0 : α1 = α2 = · · · = αP = 0.

(7)

Another way to interpret this is:
H

H0 : e ∼0 N(0, σ 2 I).

(8)

We focus on testing the significance of the first variable entering the active set and the two
main reason for doing so are as follows. First, for the following variables entering the active

14

set, the exponential distribution turns out to yield a slightly higher value than the true value
of the test statistic, so the decision tends to be conservative [47]. Second, If one would like to
include more significant regressors selected by LASSO, the Tailstop criterion [49] for ordered
selection can be used under the assumption that sequential p-values are independent. However,
we restrain from looking at p-values at greater indices because of the non-orthogonality of ξ in
our case (which leads to misleading interpretation of p-values), as pointed out in [49] and [50].
To illustrate this, we use Tailstop criteria in our simulation and it turns out that the first eighteen
regressors are selected to limit the type I error of the global null in (7) under 0.05. It is obviously
an overfitting result and Tailstop criterion cannot be trusted due to the correlation between each
row in the regressor matrix ξ . From literature, interpreting p-values from sequential hypothesis
and order selection under generic conditions is still an open question and is thus beyond the
discussion of this paper.
To test whether the first variable entering the active set is significant, we first set λ = ∞ and
gradually reduce λ until one regressor variable has a non-zero α . Denote the value of this λ as

λ1 . Also denote λ2 as the value of λ when the second regressor variable enters the active set.
For simplicity of representation, we scale the columns of the ξ so that each column has unit
norm.
Following this scaling strategy, the authors in [47] define the covariance statistic for testing
the entry of the first variable and it is written as:
T1 =

λ1 (λ1 − λ2 )
.
σ2

(9)

From Theorem 2 in [47], the authors have proved that Pr{Ti > t} → e−t as P → ∞ for all
t ≥ 0. Thus we have that the asymptotic distribution of T1 obeys exponential distribution with
parameter 1: T1 → Exp(1).
In most cases, the value of σ 2 is not known a priori. We can estimate σ 2 by the residual sum
of squared error via:

2

2

σ̂ =

α OLS )k2
k(e − ξ α̂
,
N −P

(10)

α OLS is obtained through:
where α̂
1
α OLS = arg min k(e − ξ α )k22 .
α̂
2
α

(11)

15

Plugging (11) into (9), we have a new statistic F1 , which is asymptotically following F
distribution (ratio of two independent χ 2 distribution) under the null [47]:
F1 = T1

σ 2 λ1 (λ1 − λ2 )
=
→ F2,N−P .
σ̂ 2
σ̂ 2

(12)

In conclusion, to test whether another user’s historical data can explain the residual obtained
by the univariate LASSO introduced in Section III, we compute the value of λ1 and λ2 , along
with the full linear regression in (11) which gives us σ̂ 2 . Plugging these into (12), we include
the regressor variable (which represent one particular user) entering the active set at λ1 if F1
is greater than some threshold and reject the null hypothesis in (8). Simulation results show
that performing this significance test on top of univariate LASSO type regression model has
improved the prediction by reducing relative test error by 38.3%, compared to AR(1) model.
V. S IMULATION R ESULTS
We use the data from PG&E. It contains hourly smart meter readings for residential users
during a period of one year from 2010 to 2011. Temperature data is retrieved from an online
database [51] for the same period. Some preliminary observations suggest that temperature is
not a significant regressor to predict consumption based on our data set, thus it is excluded
from the regressors for the following simulations. In addition, the consumption data has been
filtered to exclude daily periodic trends and is separated into weekday data and weekend data.
The simulation results are obtained on weekday data.
A. Prediction Evaluation Criteria
Before proceeding with the prediction methods, we first introduce the evaluation criteria for
prediction used in this paper. A naive way to evaluate prediction methods is to compare the
Mean Squared Error (MSE) or Mean Absolute Error (MAE) within the testing set:
1 n
∑ (yt − ŷt )2
n t=1

(13)

1 n
MAE = ∑ |yt − ŷt |
n t=1

(14)

MSE =

where yt is the actual value, ŷt is the predicted value, and n is the number of fitted values.

16

However, MSE and MAE are not scale invariant. This will be misleading if we want to
compare the performance for two datasets with significantly different scales and means. One
way to solve this problem is to normalize the error or to compare the relative error, i.e., the
prediction error with respect to the data scale. Here we use the MAPE mentioned earlier in
Section I, to capture the relative error:
MAPE =

1 n yt − ŷt
∑ | yt |.
n t=1

(15)

Normalized Root-Mean-Square Deviation (NRMSD) is also a good metric to compare prediction performances of data with different scales :
q
1 n
2
n ∑t=1 (yt − ŷt )
,
NRMSD =
ȳ

(16)

where ȳ is the mean of y in the training set.
In case of outliers, we adopt mean curtailing of 0.01 tail and head, or simply use the median
to replace the mean value.
In addition, to evaluate the model between good fit and complexity (in terms of the number
of parameters), Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC)
can be used [52].
B. Comparison with Different Prediction Methods
We compare our proposed method with other popular methods: 1): Simple averaging, which
averages the consumption at the same hour of the day for ten previous days; 2): LW, which
uses the consumption of the same hour of the day during last week; 3): AR(1); 4): ES with
additive error; 5): linear regression with PCA, 6): SVM with and without PCA 7): RF and 8):
NN with a single hidden layer. Input features are ten previous days’ consumption data which
has a dimension of 240, except that for NN the input features are set to be previous day’s
consumption data, which has a dimension of 24. While the first and the second method are very
intuitive methods for prediction, AR(1) and ES are popular prediction methods in analyzing
time series data. The remaining methods are popular machine learning algorithms and artificial
intelligence techniques in prediction. Particularly RF and NN are the only non linear methods
considered in this paper. They serve to provide reference results as compared to other linear
methods. Since the dataset that we use only contains time series load consumption data (plus

17

coarse zip code data), we restrain from using more complexed models which tend to overfit. In
all, our simulations using real data show that our proposed method performs similar or better to
these sophisticated algorithm while remaining a simple linear structure.
We make the following clarifications for the machine learning techniques in 5)-8). First, for
the extracted features obtained from PCA, we threshold the number of components with a
tolerance ζ on how much the chosen components explain the covariance matrix of data ξ . Then
we include those components as new features to the linear regression model. If we set up a
tolerance ζ = 0.75 and omit components if their standard deviations are smaller than ζ times
the standard deviation of the first component, we obtain an average two principle components
for each user based on our dataset. Second, a linear kernel is trained as we are comparing SVM
with several linear regression models. Third, for RF, we grow 500 trees in total. At each split of
the tree node we re-sample a third of the original features. Note that RF is a nonlinear model in
terms of input features whereas the other models are linear in input features. For a more in-depth
introduction to RF, please refer to [53]. Last, due to convergence and computational issues, we
consider one hidden layer with five neurons in the NN model.
We visualize the comparison in terms of MAPE in Fig. 2, with 50 users in the training set. In
the simulations, the forecast horizon is one hour and the shown errors are the averaged medians
over the users. Detailed comparison results based on MAPE, MAE, MSE and NRMSD are
presented in Table I, Table II, Table III and Table IV. In Table I to Table IV, the subscript W
denotes the window size. Result from linear regression after PCA is outside the range of the
results from the other prediction methods shown in Fig. 2 and is not presented. These results are
visualized in Fig 3, where a smaller error indicates a better predictive performance. In addition,
there is a constraint for the minimum number of samples in order to have a well behaved
LASSO estimator. This lower bound is of order O(slogp) [13], where s is the cardinality of the
true support and p is the total number of regressors. However, since we do not know the level of
sparsity s in advance, the number of samples should be at least O(plogp), which is around 600.
Nevertheless, since the columns in the regressor matrix are correlated, we resort to a training
size larger than the theoretical bound. In the mean time, we cannot include too many samples
into the training set to avoid losing stationarity. Here unit root test is applied to test stationarity.
To these two ends, we experiment the training sizes of 720, 960 and 1200 samples and compare
the respective performances by the proposed prediction methods.

18

0.4

MAPE

0.35

RF
NN
LASSO

ES
SVM
SVM+PCA

Averaging
LW
AR(1)

0.3

0.25

0.2
720

960

1200

Sliding window size
Fig. 2: MAPE for different prediction methods with various window size of the training set. RF
yields the best prediction performance among all aforementioned methods and LASSO yields
the best performances among linear models.

From Fig. 2, we see that RF yields the smallest prediction error in terms of MAPE. This is not
unexpected since RF is a complex method that is highly nonlinear. Surprisingly, our proposed
LASSO method matches the performance of RF when window sizes are increased. With a
sliding window with length of 720 hours, RF has an average MAPE of 0.211 and LASSO has
an average MAPE of 0.231. When we increase the size of the sliding window to include 1200
hours’ consumption data, RF has an average MAPE of 0.203 and LASSO has an average MAPE
as 0.206. At this stage, operators may prefer to use LASSO since it is as accurate as RF but has
much lower model complexity. For example, using LASSO operators can understand the most
significant hours in the past that influence the prediction. In addition, the coefficients in LASSO
can be understood as sensitivities of future prediction on past consumption, due to the linear
structure of LASSO. In contrast, the coefficients in RF do not have such “causal” structure. In

19

TABLE I: Average MAPE cross the 50 users from the dataset for the prediction methods.
Subscript W denotes the window size.
Method

MAPEW =720

MAPEW =960

MAPEW =1200

Averaging

0.258

0.258

0.258

LW

0.346

0.346

0.346

AR(1)

0.249

0.245

0.237

ES

0.250

0.244

0.236

SVM

0.327

0.282

0.264

SVM after PCA

0.305

0.304

0.285

RF

0.211

0.204

0.203

NN

0.287

0.267

0.256

LASSO

0.231

0.229

0.206

TABLE II: Average MAE cross the 50 users from the dataset for the prediction methods.
Subscript W denotes the window size.
Method

MAEW =720

MAEW =960

MAEW =1200

Averaging

0.0373

0.0373

0.0373

LW

0.0540

0.0540

0.0540

AR(1)

0.0322

0.0318

0.0311

ES

0.0340

0.0334

0.0322

SVM

0.0393

0.0338

0.0322

SVM after PCA

0.0363

0.0357

0.0368

RF

0.0273

0.0265

0.0260

NN

0.0360

0.0336

0.0311

LASSO

0.0291

0.0290

0.0260

addition, LASSO can be computed much faster than RF.
Besides these two prediction methods, SVM with either the original features or the reduced
features are average in their prediction performances. NN does not outperform other methods.

20

TABLE III: Average MSE cross the 50 users from the dataset for the prediction methods.
Subscript W denotes the window size.
Method

MSEW =720

MSEW =960

MSEW =1200

Averaging

0.00143

0.00143

0.00143

LW

0.00231

0.00231

0.00231

AR(1)

0.00103

0.00101

0.00096

ES

0.00113

0.00105

0.00094

SVM

0.00155

0.00115

0.00104

SVM after PCA

0.00132

0.00128

0.00125

RF

0.00074

0.00070

0.00067

NN

0.00123

0.00108

0.00097

LASSO

0.00085

0.00084

0.00067

TABLE IV: Average NRMSD cross the 50 users from the dataset for the prediction methods.
Subscript W denotes the window size.
Method

NRMSDW =720

NRMSDW =960

NRMSDW =1200

Averaging

0.197

0.197

0.197

LW

0.241

0.241

0.241

AR(1)

0.166

0.156

0.135

ES

0.168

0.153

0.125

SVM

0.221

0.182

0.153

SVM after PCA

0.200

0.194

0.173

RF

0.157

0.143

0.126

NN

0.178

0.167

0.165

LASSO

0.172

0.163

0.126

We believe this is due to the lack of data. Since we only have time series consumption data
without any other side information, only a shallow NN (a single hidden layer) is trained to avoid
overfitting and other computational issues.

21

0.4

0.06
0.04

MAE

MAPE

0.3
0.2
0.1

0.02

0

0

Averaging
LW
AR(1)

Averaging
LW
AR(1)

ES
SVM
SVM after PCA

RF
NN
LASSO

ES
SVM
SVM after PCA

1200
960
Window size
720

0.3

0.3

0.2

0.2

0.1

0.1

0

0

Averaging
LW
AR(1)

Averaging
LW
AR(1)

ES
SVM
SVM after PCA

RF
NN
LASSO

ES
SVM
SVM after PCA

1200
960
Window size
720

(c) Average MSE.

1200
960
Window size
720

(b) Average MAE.

NRMSD

MSE(10 -2 )

(a) Average MAPE.

RF
NN
LASSO

RF
NN
LASSO

1200
960 Window size
720

(d) Average NRMSD.

Fig. 3: Results from different error metrics for various models.

The nonparametric methods, LW and simple averaging method, return higher MAPEs compared to machine learning techniques. However, based on our dataset, simple averaging using ten
previous day’s data still outperforms LW (using data from last week only) because the noises are
leveraged by adding more consumption data into account. Since these nonparametric methods do
not depend on the size of the sliding window, the MAPEs are the same across different sliding
window sizes. What is more, these nonparametric methods are powerful when the training set
is small (a smaller sliding window size). However, when more data are available to be included
into the training set, the performances of these methods do not get improved.
In addition, the time series analysis methods, namely AR(1) and ES, yield moderate MAPE.
However, it is still not as good as either RF or LASSO. In conclusion, LASSO achieves as
good as the performance from RF, and other methods return inferior prediction performances.
Details are shown in Table I. Table II, Table III to Table IV summarize more results of the

22

LASSO
NN

Methods

RF
SVM after PCA
SVM
ES
AR(1)
LW
Averaging
0

-0.01

-0.02

-0.03

Change in MAPE
Fig. 4: Average change of MAPE with an increase of window size.

comparison based on different metrics. In all, LASSO is shown to be consistently the among
the best methods and can provide performances that are as good as advanced machine learning
techniques.
Besides the above observations on predictive performance of various models, we can also see
a difference in model sensitivity with respect to an increase in window size (a larger training
set). The results are shown in Fig 4, in terms of MAPE.
As can be seen from Fig 4, LW and averaging have no sensitivity to window size. This is
because they are both an averaged sum of a fixed history and are not parametric. AR(1) and ES
have similar sensitivities since they are both autoregressive models. SVM originally has a high
sensitivity to window size because we include many features into the model and fix a linear
kernel. However, this may not capture correctly the underlying consumption behavior. However,
with dimension reduction from PCA, the sensitivity from SVM is improved. Similarly, NN and
LASSO do not have a negligible sensitivity towards window size because of the large number

23

of features in the model but since they are more powerful models than SVM, the sensitivity is
not as high as the latter. In addition, RF falls into the scheme of ensemble learning algorithm
and its variance should be smaller than the other parametric models. Therefore its sensitivity
towards window size is small.
We then compare the averaging computational time for the aforementioned methods. The
simulations are tested on a MacBook with 2.7 GHz Intel Core i5 processor and 8GB 1867 MHz
DDR3 memory. The empirical computational time for each prediction method is shown in Table
V.
As can be seen from Table V, Simple averaging, LW, AR(1) and ES are most time efficient
but at the sacrifice of prediction accuracy. Linear regression with PCA is moderately efficient
but provides the worst prediction accuracy. RF yields the worst computational time because it
needs to be trained on the repeated bootstrapping samples. NN is also not preferable because for
a single layer model it already takes up to more than 14 seconds to train over a single sliding
window. LASSO beats both SVM with either the original regressor matrix or the principle
components in terms of time efficiency. Combining the observations from Table V with Fig. 2,
we can see that for univariate time series analysis, LASSO best trades off between prediction
accuracy and time efficiency.
Based on the observations from Fig. 2 and Table V, the complexity grows as the prediction
model varies from simple averaging method up to RF and consequently the computational time
grows accordingly. It is even more time consuming when we try to the tune the hyper-parameters
in RF, for example, the number of features to bootstrap at each split of the tree, total number
of trees, etc...What is more, although RF provides the most accuracy, it is the hardest to be
understood in terms of the estimated parameters from the model, i.e., the decision boundary at
each split of the tree node for each tree. Since throughout this paper the focus is on one hour ahead
prediction, this complexity ruins the performance of the model and is not preferred. In contrast,
LASSO yields a competing prediction performance with RF and has only one hyper-parameter
(λ ) to tune. The optimal value of this hyper-parameter can be determined from cross-validation
and is easy to implement. Thus it is the most preferable method for univariate analysis based
on our dataset.
Table VI further consolidates this preference in terms of complexity of input data. Feature
dimension is defined as the number of input variables in the model. Simple averaging method

24

has a dimension of 10 because it is a average of consumption data during the same hour for
the past ten days. LW only needs one data point during last week to make prediction. AR(1)
and ES also have a much smaller input space but their prediction performances are not as good
as either LASSO or RF. For methods associated with dimension reduction by PCA, the feature
dimension is the number of the reduced features for the model, thus around 20 in our case. We
restrict the feature dimension as 24 for NN since a higher feature dimension will require more
iterations in the algorithm to reach convergence. Other than that, for SVM, RF and LASSO,
the input data has a dimension of 240 which stands for the total number of hours in the past
ten days. With the same input complexity, LASSO yields the most efficient computational time
from Table V.
Based on the discussions from Fig. 2, Table V and Table VI, LASSO is the best choice in
terms of model efficiency and model simplicity.
TABLE V: Average computational time per user per sliding training window for each prediction
method.
Prediction method

Time (in seconds)

Averaging

less than 0.0001

LW

less than 0.0001

AR(1)

0.005

ES

0. 093

PCA + linear regression

0.175

SVM

9.347

SVM + PCA

9.522

RF

21.562

NN

14.416

LASSO

0.220

In addition, the average computational time for computing the optimal pairs using the covariance test is 3.51 seconds on average per each user. We then conclude that it is time efficient
to compute the optimal pairs of users. What is more, since the optimal pair assignment is fixed
once we compute it, this computational time can be amortized and is nearly negligible.

25

TABLE VI: Feature dimension for different prediction methods.
Prediction method

Feature dimension

Averaging

10

LW

1

AR(1)

1

ES

at most 3

PCA + linear regression

around 20

SVM

240

SVM + PCA

around 20

RF

240 (split the tree 500 times.)

NN

24

LASSO

240 (active regressors around 10-20.)

C. Results for Univariate Time Series Analysis
Based on the comparison result shown in Section V-B, in this subsection we perform more
detailed simulations on the following prediction methods from the methods considered so far:
•

Simple averaging method, which is a representative method from nonparametric methods;

•

AR(1), which is a representative method from time series methods;

•

RF, which is a representative method from machine learning techniques;

•

LASSO applied in autoregressive model, which is the proposed method in this paper.

•

SVM with and without PCA, which we use as an illustration of the failure of PCA based
on our dataset.

The simulation results associated with the covariance test are going to be introduced later
in Section V-D. The results of the prediction methods besides covariance test are shown in
Fig. 5 and Table VII, for a random selection of 150 users from the dataset. We use MAPE as
the criterion to compare prediction performances because it is a normalized metric. Again, the
MAPE shown is the median over all the sliding windows and the forecast horizon is one. As
can be seen in Fig. 5, LASSO type regression reduces the variance of MAPE for these 150
users in the dataset. Its mean is also reduced. A comparison of the aforementioned methods is

26

summarized in Table VII.
From Fig. 5, we see that MAPE varies on a user to user basis. For example, from Fig. 2(a),
using AR(1) model to predict consumption will yield a MAPE greater than 50% for some of
users. What is more, there is one extreme case for simple averaging method shown in Fig. 2(c)
that MAPE reaches 200%. This extreme value of MAPE suggests a rather unfavarable prediction
result, even worse than simply predicting zero (which has a MAPE as 100%.). However, for both
methods, a MAPE close to zero can still be achieved for other users. Besides these methods, SVM
achieves a moderate deviation of MAPE but it is still higher that by LASSO and RF. Comparing
Fig. 2(a), Fig. 2(c) with Fig. 2(b), Fig. 2(d), Fig. (e) and Fig. (f), we see that LASSO and RF
have MAPEs that are less spread between 0% and 100%. In addition, LASSO and RF have a
smaller mean of MAPE than that from the other prediction methods. More detailed statistics are
shown in Table VII. What is more, comparing Fig. 2(e) and Fig. 2(f), PCA does not improve the
performance of prediction models such as SVM, as we have discussed in Section V. B. Here,
prediction by SVM using the original features (historical consumption data) as inputs achieves
a smaller MAPE than that after dimension reduction by PCA. SVM after PCA has a MAPE of
0.291 against a MAPE of 0.264 by SVM alone.
In addition, besides the improved prediction performance compared to AR(1) and simple
averaging method, LASSO also provides some straightforward understandings of the model.
Particularly, we include as much as 240 lag orders for LASSO, which includes all the historical
data for the previous ten days. Taking user No.1 as an example, LASSO selects 16 non zero
lag orders, according to 10-fold CV with a sequence of decreasing {λk }. The lag orders that
LASSO picks are 1, 2, 5, 6, 16, 22, 23, 24, 48, 143, 144, 160, 191, 216, 238, 240. This pattern
reflects that LASSO not only selects the most recent lag orders (which is similar to AR(1)), but
the lag orders roughly at interval of one day, i.e., 24 hours as well (which is similar to simple
averaging). The coefficients for the lag orders also have different scales. In the example of user
No.1, the most recent lag orders given the largest coefficient (β1 = 0.259). The second largest
coefficient is given to lag order 24 (β24 = 0.187). The rest of the coefficients are scaled between
0.01 to 0.06.
Note that LASSO will return selected lag orders adaptively for different users. For user No.10,
the lag orders that LASSO picks up for this user are 1, 23, 24, 120, 138, 144, 168, 239, 240. In
this case, the three most recent lag orders are given the largest coefficients, where β1 = 0.579, β23

27

= 0.123 and β24 = 0.103. The rest of the coefficients scale at around 0.02. So for this particular
user, LASSO picks up less historical data at the same hour than user No. 1, at roughly one day
ago, five days ago, six days ago, seven days ago and ten days ago. This is different from simple
averaging, since LASSO does not pick up the lag order 48 or 72, which are more recent lags
than 120. Still, however different lag orders that LASSO picks for each user, we can observe a
clear pattern of multiples of 24 hours, indicating user’s consumption behavior.
In all, the proposed model based on LASSO is useful in practice because of its simplicity
and linear nature. Given the set of coefficients, utilities and planners can easily identify and
attribute the impact of different features on load consumption. For example, suppose for a given
user, our method concludes that its load at time t is mostly determined by the past load at
times t − 1, t − 3 and t − 24. This shows that this user has some short term behavior in the one
hour range, medium term behavior in the 3 hour range, and a daily cycle partner that repeats
every 24 hours. This information can be interpreted by utilities to better understand this user’s
consumption patterns and potentially identify appropriate demand programs.
D. Results for Pairing the Users
As described in Section IV, we include historical consumption data from other users to fit the
current user’s data. We compare the results of the two methods discussed as follows:
1) Reference method: LASSO selection in autoregressive model using one user’s own data.
2) Pairing by covariance test: the optimal pairing for a user is selected by the best user fitted
to the residuals obtained from this user’ consumption data after LASSO in univariate
autoregressive model. This best user corresponds to the first regressor in ξ entering the
active set when its associated p-value is small. The lag orders for this user are selected
by LASSO in the reference method, but the lag order for its paired user is fixed to one.
The results for these two methods are in Fig. 6 and Table VII, with a forecast horizon as one.
The results are based on the same 150 users chosen from the dataset. It can be seen from Fig.
6 that with the covariance test that includes another useful user into the model, the prediction
performance is enhanced. The mean of the MAPE on the test set is improved by 7% and the
variance by 6%. Note that the selection of another user into the model is not mutual, for example,
for user No.14, user No.137 is selected by covariance test whereas for user No.137, user No.50
is selected. This is essentially different from bivariate autoregressive model where pairs of users

(a). AR(1)
Frequency

100
50
0

0

0.5

1

1.5

2

MAPE
(c). Averaging
Frequency

100
50
0

0

0.5

1

1.5

2

MAPE
(e). SVM

100

Frequency

Frequency

Frequency

Frequency

28

50
0

0

0.5

1

1.5

MAPE

2

(b). LASSO

100
50
0

0

0.5

1

1.5

2

1.5

2

MAPE
(d). RF

100
50
0

0

0.5

1

MAPE
(f). SVM after PCA

100
50
0

0

0.5

1

1.5

2

MAPE

Fig. 5: Histogram comparison for univariate time series. Mean of MAPE for AR(1) is 0.339,
for LASSO 0.225, for simple averaging method 0.359, for RF 0.195 and for SVM with and
without PCA is 0.264 and 0.291 respectively. Standard deviation of MAPE for AR(1) is 0.151,
for LASSO 0.094, for simple averaging method 0.237, for RF 0.109, for SVM 0.142 and for
SVM after PCA is 0.129.

are included into the model simultaneously. Here, the pairing of the users is adaptive according
to the covariance statistic defined in (9).
As can be seen from Table VII, the proposed methods greatly improve the prediction performance by reducing the relative prediction error, thus will provide more reliable results on
consumption prediction. The proposed prediction methods in this paper can help compare the
user behavior with and without demand response incentives. In the future, we will use the

29

TABLE VII: Results for the prediction methods, where CovT stands for covariance test and Ave
stands for simple averaging method.
Ave

AR(1)

LASSO

RF

SVM

CovT

mean (MAPE)

0.359

0.339

0.225

0.195

0.264

0.209

sd (MAPE)

0.237

0.151

0.094

0.109

0.142

0.088

Frequency

(a). LASSO
40
20
0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.4

0.5

0.6

Frequency

MAPE
(b). CovT
50

0

0

0.1

0.2

0.3

MAPE
Fig. 6: Histogram comparison for multivariate time series, where covT stands for covariance test.
Mean MAPE for LASSO is 0.225 and for performing covariance test after LASSO is 0.209.
Standard deviation of MAPE has also reduced by 7% from LASSO to performing covariance
test after LASSO.

results in this paper to investigate the demand response effect of users by comparing their actual
consumption data and predicted consumption data.
VI. C ONCLUSIONS

AND

F UTURE WORKS

In this paper, we propose a sparse autoregressive model to predict the electricity consumption
of individual users. First, we adopt LASSO to recover the sparsity in linear regression models.
LASSO selects the most recent lag orders and important lag orders with close multiples of
24 hours, which reveals user consumption patterns. Second, we improve prediction accuracy

30

of a particular user by leveraging other user’s historical data, where we use the covariance
statistic to test if inclusion of another user is significant to explain the fitted residual from
univariate LASSO-type regression. Third, we give a comprehensive analysis on the proposed
methods against several available predictive models in literature. We observe that LASSO best
trades off between model complexity and prediction performance. It has the least predictive
error among the linear models and yields a competitive performance as compared to ensemble
methods such as random forests. Further simulation results show that LASSO with covariance
test outperforms both simple averaging method and auto regressive method with order 1 for
individual consumption prediction, improving the prediction performance by 38.4% in terms of
MAPE.
R EFERENCES
[1] H. Afares and M. Nazeeruddin, “Electric load forecasting: literature survey and classification of methods,” International
Journal of Systems Science, vol. 33, no. 1, pp. 23–34, 2012.
[2] J. W. Taylor and P. E. McSharry, “Short-term load forecasting methods: An evaluation based on european data,” Power
Systems, IEEE Transactions on, vol. 22, no. 4, pp. 2213–2219, 2007.
[3] R. Weron, Modeling and forecasting electricity loads and prices: a statistical approach. John Wiley & Sons, 2007, vol.
403.
[4] R. Sevlian and R. Rajagopal, “Short term electricity load forecasting on varying levels of aggregation,” arXiv preprint
arXiv:1404.0058, 2014.
[5] M. Baran and T. McDermott, “Distribution system state estimation using ami data,” in 2009 Power Systems Conference
and Exposition.

IEEE, 2009, pp. 1–3.

[6] C. Chen, J. Wang, Y. Heo, and S. Kishore, “MPC-based appliance sheduling for residential building energy management
controller,” Smart grid, IEEE Transactions on, vol. 4, no. 3, pp. 1401–1410, 2013.
[7] M. F. S. Aman, C. Charalampos, U. Noor, Y. Simmhan, and V. Prasanna, “Empirical comparison of prediction methods
for electricity consumption forecasting,” University of Southern California, Tech. Rep., 2014.
[8] J. MacDonald, P. Cappers, D. Callaway, and S. Kiliccote, “Demand response providing ancillary services,” Grid-Interop,
Tech. Rep., 2012.
[9] P. Li, X. Guan, J. Wu, and D. Wang, “An integrated energy exchange scheduling and pricing strategy for multi-microgrid
system,” in In TENCON 2013-2013 IEEE Region 10 Conference (31194).

IEEE, 2013, pp. 1–5.

[10] D. Wang, X. Guan, J. Wu, P. Li, P. Zan, and H. Xu, “Integrated energy exchange scheduling for multimicorgrid system
with electric vehicles,” Smart grid, IEEE Transactions on, vol. 7, no. 4, pp. 1762 – 1774, 2016.
[11] R. Tibshirani, “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society. Series B
(Methodological), pp. 267–288, 1996.
[12] P. Zhao and B. Yu, “On model selection consistency of lasso,” The Journal of Machine Learning Research, vol. 7, pp.
2541–2563, 2006.

31

[13] M. J. Wainwright, “Sharp thresholds for high-dimensional and noisy sparsity recovery using-constrained quadratic
programming (lasso),” Information Theory, IEEE Transactions on, vol. 55, no. 5, pp. 2183–2202, 2009.
[14] H. Zou, “The adaptive lasso and its oracle properties,” Journal of the American statistical association, vol. 101, no. 476,
pp. 1418–1829, 2006.
[15] L. Meier, S. V. D. Geer, and P. Bhlmann, “The group lasso for logistic regression,” Journal of the Royal Statistical Society:
Series B (Statistical Methodology), vol. 70, no. 1, pp. 53–71, 2008.
[16] H. Xu, D. Eis, and P. Ramadge, “The generalized lasso is reducible to a subspace constrained lasso,” in 2013 IEEE
International Conference on Acoustics, Speech and Signal Processing.

IEEE, 2013, pp. 3268–3272.

[17] F. Bunea, J. Lederer, and S. Y., “The group square-root lasso: Theoretical properties and fast algorithms,” Information
Theory, IEEE Transactions on, vol. 60, no. 2, pp. 1313–1325, 2014.
[18] L. Breiman, “Random forests,” Machine learning, vol. 45, no. 1, pp. 5–32, 2001.
[19] T. K. Ho, “Random decision forests,” in Document Analysis and Recognition, 1995 Proceedings of the Third International
Conference on, vol. 1. IEEE, 1995, pp. 278–282.
[20] F. Ziel, R. Steinert, and S. Husmann, “Efficient modeling and forecasting of electricity spot prices,” Energy Economics,
vol. 47, pp. 98–111, 2015.
[21] T. Hong, “Energy forecasting: Past, present, and future,” Foresight: The international Journal of Applied Forecasting,
vol. 32, pp. 43–48, 2014.
[22] NewYorkISO, “Caiso demand response resource user guide,” NewYorkISO, Tech. Rep., 2007.
[23] C. Katie, M. A. Piette, C. A. Goldman, and S. Kiliccote, “Statistical analysis of baseline load models for non-residential
loads,” Energy and Buildings, vol. 41, no. 4, pp. 374–381, 2009.
[24] S. D. Ramchurn, P. Vytelingum, A. Rogers, and N. R. Jennings, “Putting the smarts into the smart grid: a grand challenge
for artificial intelligence,” Communications of the ACM, vol. 55, no. 4, pp. 86–97, 2012.
[25] S. Aman, Y. Simmhan, and V. K. Prasanna, “Improving energy use forecast for campus microgrids using indirect indicators,”
in Data Mining Workshops (ICDMW), 2011 IEEE 11th International Conference on.

IEEE, 2011, pp. 389–397.

[26] B. J. Chen, M. W. Chang, and C. J. Lin, “Load forecasting using support vector machines: a study on eunite competition
2001,” Power Systems, IEEE Transactions on, vol. 19, no. 4, pp. 1821–1830, 2004.
[27] Y. Weng and R. Rajagopal, “Probabilistic baseline estimation via gaussian process,” in 2015 IEEE Power and Energy
Society General Meeting.

IEEE, 2015, pp. 1–5.

[28] H. Quan, D. Srinivasan, and A. Khosravi, “Short-term load and wind power forecasting using nueral network-based
prediction intervals,” Neural Networks and Learning Systems, IEEE Transactions on, vol. 25, no. 2, pp. 303–315, 2014.
[29] A. Sarwat, M. Amini, A. Domijan, A. Damnjanovic, and F. Kaleem, “Weather-based interruption prediction in the smart
grid utilizing chronological data,” Journal of Modern Power Systems and Clean Energy, vol. 4, no. 2, pp. 308–315, 2016.
[30] J. B. Fiot and F. Dinuzzo, “Electricity demand forecasting by multi-task learning,” Smart Grid, IEEE Transactions on,
vol. PP, no. 99, pp. 1–8, 2016.
[31] O. Abedinia, N. Amjady, and H. Zareipour, “A new feature selection technique for load and price forecast of electrical
power systems,” Power Systems, IEEE Transactions on, vol. PP, no. 99, p. 1, 2016.
[32] N. Amjady, “Short-term hourly load forecasting using time-series modeling with peak load estimation capability,” Power
Systems, IEEE Transactions on, vol. 16, no. 3, pp. 498–505, 2001.
[33] S. Chan, K. M. Tsui, H. C. Wu, Y. Hou, Y. C. Wu, and F. Wu, “Load/price forecasting and managing demand response
for smart grid: Methodologies and challenges,” Signal Processing Magazine, IEEE, vol. 29, no. 5, pp. 68–85, 2012.

32

[34] J. W. Taylor, “Short-term load forecasting with exponentially weighted methods,” Power Systems, IEEE Transactions on,
vol. 27, no. 1, pp. 458–464, 2012.
[35] M. Hagan and S. Behr, “The time series approach to short term load forecasting,” Power Systems, IEEE Transactions on,
vol. 2, no. 3, pp. 785–791, 1987.
[36] M. H. Amini, O. Karabasoglu, M. D. Ilic, K. Boroojeni, and S. S. Iyengar, “Arima-based demand forecasting method
considering probabilistic model of electric vehicles’ parking lots,” in 2016 IEEE Power and Energy Society General
Meeting.

IEEE, 2015, pp. 1–5.

[37] J. Contreras, R. Espinola, F. Nogales, and A. Conejo, “Arima models to predict next-day electricity prices,” Smart Grid,
IEEE Transactions on, vol. 18, no. 3, pp. 1014–1020, 2003.
[38] M. Amini, A. Kargarian, and O. Karabasoglu, “Arima-based decoupled time series forecasting of electric vehicle charging
demand for stochastic power system operation,” Electric Power Systems Research, pp. 378–390, 2016.
[39] G. Chicco, “Overview and performance assessment of the clustering methods for electrical load pattern grouping,” Energy,
vol. 42, no. 1, pp. 68–80, 2012.
[40] J. D. Hamilton, Time series analysis.

Princeton university press, 1994.

[41] J. Friedman, T. Hastie, and R. Tibshirani, The elements of statistical learning.

Springer, 2001.

[42] N. Rao, R. Nowak, C. Cox, and T. Rogers, “Classification with the sparse group lasso,” Signal Processing, IEEE
Transactions on, vol. 64, no. 2, pp. 448–463, 2016.
[43] A. Jung, G. Hannak, and N. Goertz, “Graphical lasso based model selection for time series,” Signal Processing Letters,
IEEE, vol. 22, no. 10, pp. 1781–1785, 2015.
[44] D. Angelosante, J. A. Bazerque, and G. B. Giannakis, “Online adaptive estimation of sparse signals: Where rls meets
the-norm,” Signal Processing, IEEE Transactions on, vol. 58, no. 7, pp. 3436–3447, 2010.
[45] N. Ludwig, S. Feuerriegel, and D. Neumann, “Putting big data analytics to work: Feature selection for forecasting electricity
prices using the lasso and random forests,” Journal of Decision Systems, vol. 24, no. 1, pp. 19–36, 2015.
[46] Y. Goude, R. Nedellec, and N. Kong, “Local short and middle term electricity load forecasting with semi-parametric
additive models,” Smart Grid, IEEE Transactions on, vol. 5, no. 1, pp. 440–446, 2014.
[47] R. Lockhart, J. Taylor, R. J. Tibshirani, and R. Tibshirani, “A significance test for the lasso,” Annals of statistics, vol. 42,
no. 2, p. 413, 2014.
[48] R. J. Tibshirani, “The lasso problem and uniqueness,” Electronic Journal of Statistics, vol. 7, pp. 1456–1490, 2013.
[49] M. G. G’Sell, S. Wager, A. Chouldechova, and R. Tibshirani, “Sequential selection procedures and false discovery rate
control,” arXiv preprint arXiv:1309.5352, 2013.
[50] A. Buja and L. Brown, “Discussion: a significance test for the lasso,” The Annals of Statistics, vol. 42, no. 2, pp. 509–517,
2014.
[51] National Oceanic and Atmospheric Administration, “National Centers for Environmental Information,” 2015. [Online].
Available: http://cdo.ncdc.noaa.gov
[52] K. G. Boroojeni, M. H. Amini, S. Bahrami, S. S. Iyengar, A. I. Sarwat, and O.Karabasoglu, “A novel multi-timescale modeling for electric power demand forecasting: from short-term to medium-term horizon,” Electric Power Systems
Research, vol. 142, pp. 58–73, 2017.
[53] A. Liaw and M. Wiener, “Classification and regression by randomforest,” R news, vol. 2, no. 3, pp. 18–22, 2002.

33

A PPENDIX A
S TATISTICS

CONCEPTS

Since many statistical terms appear throughout the paper, we therefore give more explanations
to them in this appendix.
A. Univariate autoregressive model
An univariate autoregressive model of consumption data yt is expressed in the following form:
I

yt = β0 + ∑ βi yt−i + εt .

(17)

i=1

where yt is a scalar at each time t. It simply refers to the autoregressive models for one
dimensional time series data.
B. Vector autoregressive model
A vector autoregressive model aims to jointly predict multiple users’ future demand by their
past consumption data. It is expressed in the following form:
I

zt = A 0 + ∑ A i z t−i + ε t .

(18)

i=1

where zt is a consumption vector at each time t for multiple users. Ai is the parameter matrix at each lag order i associated with these users’ past consumption data. Therefore, vector
autoregressive model refers to the autoregressive model for multi-dimensional time series data.
C. Out-of-sample error
Out-of-sample error depicts prediction error of the statistical model. More specifically, the
out-of-sample error is the deviation from consumption prediction to the real consumption that
is not used to train the model by the operator. For example, in OLS estimation, the objective is
to minimize the sum of squared error for the samples in the training set. Out-of-sample error is
thus the estimation error in the test set.

34

D. AIC/BIC
AIC (Akaike information criterion) is a measure of relative performance of statistical models.
It is similar to traditional error metrics such as MSE only that model complexity is involved in
its computation. AIC computed as:
AIC = 2k − ln L

(19)

where k is the number of free parameters in the model and L is the likelihood.
Similarly, BIC (Bayesian Information Criterion) is also a measure of relative performance of
statistical models. It is defined as the following:
BIC = −2 ln L + k ln n

(20)

where n is the number of data points in the dataset.
E. Active set
An active set is the set of all non zero coefficients estimated at a fixed value of the tuning
parameter λ in LASSO. In this paper, the active set for one particular user includes the indices
of users that are potentially useful in predicting this user’s future demand.
F. Tailstop criterion
Tailstop criterion is introduced in [49] and is used as a stopping criterion in a sequence of
significance tests. For more details, please refer to [49].
G. Covariance statistic
Covariance statistic is used in the significance test for LASSO. More specifically, it is the
statistic that tests whether the inclusion of the variables are significant or not. In this paper, it
is used to select the index of the user whose past consumption data is most useful in predicting
one particular user’s future demand. The larger the statistic, the more significant the selected
user’s influence is on the prediction.

35

A PPENDIX B
M ODEL

COMPARISON USING

AIC/BIC

Besides the error metrics introduced in Section V. A, AIC and BIC are two other criteria
to model selection. We therefore use AIC and BIC to compare the autoregressive models with
order up to 5 and the proposed sparse model obtained by LASSO. Since the objective for both
models is to minimize the RSS (Residual Sum of Squares), AIC/BIC can be reduced to:
AIC = 2k + n ln(RSS) +C

(21)

where C is a constant across models, and
BIC = n ln(RSS/n) + k lnn

(22)

For autoregressive models, k is equal to the number of lag orders in the model. For LASSO,
k is defined as the level of sparsity of the model.
A comparison between AR(x) model and LASSO based on AIC/BIC is provided in Table VIII
and in Fig. 7, with a randomly selected user from the pool. Note that lower values of AIC and
BIC are preferred, however their absolute value does not indicate whether the model is good or
not. Therefore AIC and BIC criteria are relative measures of the model quality.

3650

-4850

3600

-4900

3550

-4950

3500

-5000

3450

-5050

3400

-5100

3350
3300
3250
AR(1)

-5150
BIC
Relative AIC

AR(2)

-5200
AR(3)

AR(4)

AR(5)

-5250
LASSO

Fig. 7: Comparison between different autoregressive models using AIC and BIC.

36

TABLE VIII: Comparison between different autoregressive models using AIC and BIC.
Method

AIC

BIC

AR(1)

3645.8+C

-4892.7

AR(2)

3645.0+C

-4888.4

AR(3)

3594.4+C

-4933.8

AR(4)

3575.1+C

-4950.0

AR(5)

3560.1+C

-4957.7

LASSO

3283.5+C

-5203.0

As can be seen from Table VIII and Fig. 7, LASSO has both the lowest AIC as well as BIC,
which implies that LASSO is the best model among the models presented in Table VIII. From
Table VIII we observe that the AIC and BIC from AR(x) models are quite similiar, with only a
change of 1.3%. However, AR(5) has an increase of 23.8% in out-of-sample error, as shown in
Fig. 1. Therefore AIC/ BIC is not a preferable criteria to compare predictive performance based
on our data set.

