Wavelet Estimation for Factor Models with Time-Varying Loadings

by   Duván Humberto Cataño, et al.

We introduce a high-dimensional factor model with time-varying loadings. We cover both stationary and nonstationary factors to increase the possibilities of applications. We propose an estimation procedure based on two stages. First, we estimate common factors by principal components. In the second step, considering the estimated factors as observed, the time-varying loadings are estimated by an iterative generalized least squares procedure using wavelet functions. We investigate the finite sample features by some Monte Carlo simulations. Finally, we apply the model to study the Nord Pool power market's electricity prices and loads.



page 19


Modelling Time-Varying First and Second-Order Structure of Time Series via Wavelets and Differencing

Most time series observed in practice exhibit time-varying trend (first-...

Joint Mean-Vector and Var-Matrix estimation for Locally Stationary VAR(1) processes

During the last two decades, locally stationary processes have been wide...

Time varying regression with hidden linear dynamics

We revisit a model for time-varying linear regression that assumes the u...

Sparse time-varying parameter VECMs with an application to modeling electricity prices

In this paper we propose a time-varying parameter (TVP) vector error cor...

Inferential Theory for Granular Instrumental Variables in High Dimensions

The Granular Instrumental Variables (GIV) methodology exploits panels wi...

Determining the dimension of factor structures in non-stationary large datasets

We propose a procedure to determine the dimension of the common factor s...

Uniform Inference for Conditional Factor Models with Instrumental and Idiosyncratic Betas

It has been well known in financial economics that factor betas depend o...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Factor models have been widely used in the last decade due to their ability to explain the structure of a common variability among time series by a small number of unobservable common factors. In this sense, these models are used to reduce the dimensionality of complex systems. The studies of factor models have encompassed the stationary and nonstationary frameworks by different estimation methods, see among many others, pena1987identifying, forni2000generalized, stock2002forecasting, bai2002determining, bai2003inferential, forni2004generalized, forni2005generalized, lam2012factor, bai2016maximum, chan2017factor, and gao2019structural for stationary cases, and bai2004panic, pena2006nonstationary, barigozzi2016non, and rodriguez2017estimation for nonstationary cases.

The literature has recently focused on a particular framework of factor models that allow the loadings to vary over time. motta2011locally introduce deterministic smooth variations in factor loadings and propose estimation procedures based on locally weighted generalized least squares using kernel functions under a time-domain approach. Similarly, eichler2011fitting allow for a nonstationary structure in the factor model with deterministic time-dependent functions in factor loadings. Their estimation procedure can be seen as a time-varying spectral density matrix of the underlying process. Furthermore, mikkelsen2018consistent

propose a factor model with time-varying loadings that evolve as stationary VAR processes employing Kalman filter procedures to obtain the maximum likelihood estimators of the factor loadings parameters.

In this paper, we use wavelet functions to define smooth variations of loadings in high-dimensional factor models. Our model can be helpful in some economic applications when the dynamics are driven by smooth variations, whose cumulative effects cannot be ignored, see suwang2017, and baihan2016

, for details. In our framework, to capture the common smooth variations in the vector of time series, parameters in the loading matrix are assumed to be well approximated by deterministic functions over time.

The estimation procedure consists of two steps. First, following mikkelsen2018consistent

, common factors are estimated by principal component analysis (PCA). In the second step, using estimated factors of the first stage, factor loadings are estimated by an iterative procedure that combines generalized least squares (GLS) using wavelet functions. We show that factors estimated by principal components are consistent after controlling the magnitude of the loadings’ instabilities. We highlight that a requirement for such consistency is that common factors need to be independent of factor loadings’ functions. We also study the case of nonstationary factors following the approach of

pena2006nonstationary, who propose to use a generalized covariance matrix to estimate common factors.

We use Monte Carlo simulations to show that the method correctly identifies the factors and loadings even in relatively small samples. Finally, we use the methodology proposed to analyze the Nord Pool power market. We use the electricity system prices and loads throughout two years and a half to show that factor loadings are not invariant along time, illustrating our model’s usefulness. We find that some features of electricity prices and loads, see, e.g., weron2007modeling, and Weron20141030, are well extracted by common factors estimates and by time-varying loadings estimates.

In short, the contributions of this article are i) We propose wavelet functions to estimate time-varying loadings and a consistent method to estimate them; ii) we also allow for nonstationary factors to cover more possibilities of application, and iii) the usefulness of the model is discussed through a Monte Carlo study, and by the empirical application.

The remainder of the paper is organized as follows. Section 2 introduces the model, shows the consistency of principal components with stationary and nonstationary variables, and introduces the wavelet functions. Section 3 discusses the estimation procedure. Monte Carlo experiments are presented in Section 4, whereas an empirical illustration is provided in Section 5. Finally, Section 6 concludes. We use bold-unslanted letters for matrices, bold-slanted letters for vectors, and unbold letters for scalars. We denote by the trace operator, by the rank of a matrix , by

the identity matrix of dimension

, by the Kronecker product and by the Frobenius (Euclidean) norm, i.e.,

2 The model

In this section, we introduce the model. First, we consider stationary factors to motivate the general setup and the estimation procedure. Then, in the next section, we relax the stationarity assumption to allow for nonstationary factors.

2.1 Stationary factors

The model we propose is as follows


where the common component, , is an dimensional locally stationary process, in the sense of dahlhaus1997fitting, and the loadings are defined in the re-scaled time, . are the unobservable common factors and , the idiosyncratic component, is a sequence of weakly dependent variables, see e.g. stock2002forecasting, and bai2003inferential.

, is the time-dependent factor loading matrix. In this respect, the influence of on the observed process varies over time, and captures the smooth variations of the loadings with respect , a constant loading matrix.

When , the standard factor model is considered; therefore, the model proposed in (1) can be seen as a generalization of the standard factor model of bai2003inferential.

Definition 1.

The sequence in (1) follows a factor model with time-varying loadings if:

  1. For , there is a function with

    such that

    with and is a positive definite diagonal matrix.

  2. is a positive definite matrix.

With an arbitrary constant , we then use the following assumptions to study the model in (1) following bai2002determining.

Assumption A.


  1. [label=]

  2. as with is a positive definite diagonal matrix.

Assumption B.

Factor loadings:

  1. [label=]

  2. and , as , with positive definite matrix , and is the th row of

  3. where is the th row of

  4. ,  for and

Assumption C.

Idiosyncratic terms:

  1. [label=]

  2. with , for a constant and

  3. and

Assumption D.

Time-varying factor loadings and factors:
, and are functions such that the following conditions are fulfilled for any and

  1. [label=]

  2. ,

  3. .

Assumption E.

The process and are independent of each other for any .

Assumption A imposes standard moment conditions. The unobservable factors have finite fourth moments, and their covariance converges in probability to a positive definite matrix. Assumptions


ensure that each factor has a nontrivial contribution to the variance of

. Assumption ensures the existence of the expansion in the wavelet function for the factor loadings. Assumption C allows dependence in the idiosyncratic process. Assumption D is required to guarantee the consistency of the principal components. Finally, independence between factors and the idiosyncratic term is provided in Assumption E.

2.1.1 Principal component estimators

We use PCA to estimate the factors . It is well-known that principal components are obtained solving the optimization problem


where is a matrix and is a matrix. We need to impose some restrictions to guarantee the identification of the parameters. Solving for , the normalization provides the necessary number of restrictions.

With this, minimize (3) is equivalent to maximize , where . Then, the estimated factor matrix, , is

times the eigenvectors corresponding to the

largest eigenvalues of the

matrix . It is well-known that this solution is not unique, that is, any orthogonal rotation of is also a solution. See bai2008large for more details.

The following theorem is a modified version of Theorem 1 in bates2013consistent and shows that, under the assumptions previously stated, it is possible to consistently estimate any rotation of the factors by principal components even if the loadings are time-varying.

Theorem 1.

Under Assumptions A-E, there exists an matrix , such that


as , where with , , and defined in the Assumption D. Furthermore,

where is a diagonal matrix of the largest eigenvalues of the matrix .

Proof. See the appendix A.1.

Theorem 1 points out that the average squared deviation between the estimated factors and the space spanned by a rotation of the actual factors will vanish at rate , which is similar to that in bai2002determining. Note that from (4), the estimated common factors, , are identified through a rotation, then, principal components converge to a rotation of the actual common factors .

2.2 The nonstationary model

In many areas, as in economics and finances, finding strong evidence of nonstationary processes has repeatedly been reported in many empirical studies. Consequently, it is natural to think that many panel data may include nonstationary economic or financial variables. There has been some debate concerning the use of differenced variables. The main argument being discussed is whether differencing the series causes a severe loss of information.

In this section, we allow for nonstationary variables in model 1. Then, assume that is . Note that while assumptions B-E do not differ under the nonstationary setup, we should modify the factor structure’s assumption. With this in mind, we follow the assumption 1 in pena2006nonstationary to define the nonstationary common factors as follows

Assumption F.

Nonstationary Factors:

where is the lag operator, is a positive integer, is a vector of drifts, and for any matrix or vector M. Define with

To estimate the nonstationary model, we use the methodology proposed by pena2006nonstationary, who use generalized covariance matrices defined as


where and can be either 0 or 1.

Under this framework, consistent factor estimates are given by


where is a matrix, is a matrix and is a matrix composed by the first eigenvectors of . Following the same reasoning as in Theorem 1, it can be shown that the average squared deviation between the estimated factors (6) and space spanned by a rotation of the actual factors will vanish as

2.3 Wavelets

The basic idea of a wavelet is to construct infinite collections of translated and scaled versions of the scaling function and the wavelet such as , and for . Suppose that forms an orthonormal basis of , for any coarse scale . A key point is to construct and with a compact support that generates an orthonormal system, which has location in time-frequency. From this, we can get parsimonious representations for a wide class of wavelet functions, further details in changmore2005, and portomoreaub2008.

In some applications, these functions are defined in a compact set in for functions , for and defined in (2). Then, it is necessary to consider an orthonormal system that generates . For the construction of these orthonormal systems, we follow the procedure by cohen1995wavelets, that generates multiresolution levels , where the spaces are generated by Negative values of are not necessary since , and if , , see vidakovic2009statistical for more details and different approaches. Therefore, for any function , can be expand in series of orthogonal functions


where we take and . For each , the set generates values of such that belongs to the scale . For example, for , there are eight wavelet coefficients in the scale , whereas for , only four coefficients in the scale .

Some applications consider the equation in (7) for a maximum resolution level , through


In this way, the function approximates to the space . We use ordinary wavelets as in dahlhaus1997identification because the performance is suitable in the case of smooth functions. Particularly, we employ Daubechies (, hereafter) and Haar wavelets of compact supports.

3 Estimation of time-varying loadings by wavelets

We consider the process in (1) with common factors to discuss the estimation procedure of the time-varying loadings. From now on, we consider the loading matrix, , as a unique function over time , given by


with , and In matrix form, we have


where is an dimensional vector of time series, is the time-varying loading matrix with , for , and . is the common factor, and is the idisyncratic process. From (9), and the Assumption E, the structure of the covariance matrix of the process is written as

that means, the variance of the common component is .

For the construction of the time-varying loadings estimates, we first assume that the estimator of the common factors are obtained by PCA of the dimensional time series . Then, functions of time-varying loadings are approximated in series of orthogonal wavelets as in (8), for a fixed resolution level ,


The values of vary depending on the resolution level in the wavelet decomposition. We choose the maximum resolution , such that , see dahlhaus1997identification for details of this selection. In practice, the coefficients are obtained for a particular estimation method. In this paper, we use GLS to estimate these coefficients and to reconstruct the loadings functions.

Let be time series with which are generated by


where the common factors, , are estimated by principal components. Each loading function is written as in (11), then when plugging each into (10), we have



are matrices for and is null matrix.

Let be a matrix, then

is matrix that depends on the estimated factors, , the wavelets , and the resolution level , with vector of parameters of dimension for , where .

Each is composed by the wavelets coefficients of the th row of the matrix . Therefore, the total number of wavelets parameters to be estimated is .

Hence, the model in (13) can be represented in a linear model form as

where is the response vector and

is the usual design matrix in regression analysis. Assuming that the covariance matrix of the idiosyncratic errors,

, is known, then the GLS estimator of the coefficients is given by

where is a matrix defined as

These results provide linear estimators of wavelet coefficients for the time-varying loadings assuming that the covariance matrix of the idiosyncratic error is known. Some procedures, such as maximum likelihood methods, are not computationally efficient for estimating such a model since the number of parameters tends to be very large. In this light, we use GLS to simplify the implementation.

Note that we can use a different basis of wavelet functions , and for each . In the simulation section, we use a similar basis to simplify the exposition.

3.1 Estimation algorithm

The estimation algorithm is as follows:

Step 1.

Use principal components to estimate as

where is the eigenvector corresponding to the th largest eigenvalue, for , of the matrix . Here, and denote and matrices, respectively. Next, the loading function matrix, , is approximated by wavelets

with , , and . Then, we write the equation (12) as

where and are vectors, and the dimensions of and are and , respectively, where indicates the resolution level chosen in the wavelet expansions.

Step 2.

Estimate by GLS the wavelet coefficients as

where , , and is used as initial value.

Step 3.

Using the estimated coefficient in the Step 2, the loadings are obtained as


Step 4.

With , obtain residuals, Then, compute

Step 5.

Back to step 2 with . Iterate -times the procedure to obtain the sequences . Stop the iteration when

for any small , where denotes the Frobenius norm for

It is possible to improve the efficiency of the common factors by giving the new value of the loading matrix, regressing the series on loadings to obtain a new estimate of the factors, and then iterate the algorithm. However, this procedure is computationally very intensive for large data sets. As shown in the following two sections, it is unnecessary to implement an expensive iterative procedure to get successful results. Our algorithm is comparable with the methodology proposed by mikkelsen2018consistent, which also does not iterate the algorithm.

4 Monte Carlo simulation

We examine the finite-sample properties of the estimation procedure proposed above using a Monte Carlo study. The model in (1) is generated as

where is the lag operator, , for , and the matrix is generated by two structures; i) , a Toeplitz matrix, and ii) a diagonal matrix. Furthermore, at time , denotes the th time series, is vector of loadings, which are alternately generated by some smooth functions. We discuss a couple of functions used below. is the vector of factors, and factor errors and are vectors of idiosyncratic terms which are independent to each other.

A referee kind let us know that Gaussian distributions for both error terms in our Monte Carlo experiment can be relaxed by incorporating other types of distributions. The referee points out that assuming homogeneous distributions on

can help the model be more accurate in some empirical applications. As commented by the referee, a possible estimation method can be the wavelet polynomial chaos method.

Another possibility to explore our model under a non-Gaussian distribution could be under the state-space modeling. As pointed in durbin2012time

, it is common to assume normality distribution in the innovation in state-space models because the model is estimated by maximizing a Gaussian log-likelihood, which is evaluated by the Kalman filter. In principle,quasi-maximum likelihood methods can be used when actual distributions of the error terms are non-Gaussian. The possibility of non-Gaussian distributions is beyond the scope of the present paper and is not further explored.

poncela2021factor provide an excellent survey on factor extraction using Kalman filter.

In our Monte Carlo study, the model is generated with cross-sectional units and sample sizes. We consider for simplicity only two common factors . Furthermore, three values for for , are considered. Two values for ; a Toeplitz matrix with for correlated noise and

for uncorrelated, where a uniform distribution

generates the entries of the diagonal matrix. Note that simulated data are standardized before extracting the principal components. Common factors are estimated by principal components in cases with , and by the procedure of pena2006nonstationary in the case with . All simulations are based on 1000 replications of the model.

We rotate the obtained factors to compare proposed estimations with the actual factors. The optimal rotation is obtained by maximizing
. The solution is given by where and are orthogonal matrices of the decomposition . When the number of principal components is not equal to the number of factors , we rotate the first principal components, see eickmeier2015classical

. Both estimated and simulated factors are re-scaled to keep the same standard deviation, then


where is the th column of the matrix of the rotated principal components .

As explained before, these rotated factors are now treated as observed variables in the regression model

to estimate the wavelet coefficients, , where .

To investigate the performance of the estimation procedure, estimated and simulated factors are compared as follows:

  1. The precision of the estimation factors is measured by the statistics as in bates2013consistent, given by


    where is the matrix of estimated factors as in (14) and is the matrix of the actual factors, that is the simulated ones. This statistics is a multivariate in a regression of the actual factors on the principal components. When the canonical correlation of the actual and estimated factors tends to one, then as well.

  2. We measure the precision of loadings estimates by mean square errors (MSE) between estimated and actual loadings, as in motta2011locally. The MSE is computed as follows

    for . The estimator of the factor loadings matrix, , is chosen by a path such that

Wavelet functions Haar D8 20 512 0.9341 0.0103 0.8042 0.1970 0.7450 0.1562 0.9047 0.0111 0.8655 0.0023 0.7314 0.1191 1024 0.9421 0.0076 0.8135 0.0900 0.7501 0.1202 0.9347 0.0082 0.8764 0.0011 0.7439 0.0876 2048 0.9432 0.0018 0.8237 0.0052 0.7840 0.0711 0.9470 0.0005 0.8855 0.0009 0.8046 0.0045 30 512 0.9236 0.0098 0.8056 0.0128 0.7091 0.1091 0.9460 0.0071 0.8521 0.0127 0.7931 0.1093 1024 0.9547 0.0051 0.8125 0.0090 0.7112 0.0859 0.9487 0.0047 0.8671 0.0079 0.7963 0.0759 2048 0.9723 0.0015 0.8268 0.0027 0.7324 0.0738 0.9498 0.0020 0.8691 0.0059 0.8213 0.0028 100 512 0.9062 0.0088 0.8019 0.0147 0.7014 0.1198 0.9643 0.0062 0.8674 0.0102 0.8984 0.0993 1024 0.9100 0.0042 0.8418 0.0081 0.7020 0.0781 0.9701 0.0037 0.9044 0.0081 0.9105 0.0651 2048 0.9241 0.0017 0.8857 0.0035 0.7153 0.0719 0.9812 0.0025 0.9202 0.0054 0.9300 0.0198 20 512 0.8631 0.0672 0.6932 0.0891 0.6711 0.2201 0.8911 0.0128 0.6871 0.0593 0.6242 0.1223 1024 0.8634 0.0501 0.6953 0.0702 0.7012 0.1901 0.8926 0.0100 0.7001 0.0337 0.6832 0.0693 2048 0.8911 0.0137 0.7723 0.0433 0.7321 0.1642 0.8971 0.0091 0.7307 0.0108 0.6941 0.0031 30 512 0.8531 0.0472 0.7984 0.0621 0.6812 0.0912 0.9021 0.0117 0.6815 0.0311 0.7483 0.1114 1024 0.8671 0.0231 0.7730 0.0539 0.7122 0.0734 0.9126 0.0090 0.7270 0.0276 0.7126 0.0554 2048 0.8711 0.0092 0.8200 0.0311 0.7212 0.0819 0.9232 0.0068 0.7305 0.0109 0.7531 0.0037 100 512 0.8503 0.0398 0.8081 0.0510 0.7294 0.0831 0.9368 0.0101 0.7902 0.0263 0.7566 0.0911 1024 0.8602 0.0288 0.7942 0.0495 0.7413 0.0698 0.9404 0.0087 0.8012 0.0201 0.7774 0.0283 2048 0.8893 0.0104 0.8491 0.0209 0.7629 0.0503 0.9499 0.0070 0.8125 0.0117 0.7849 0.0325

Notes: The DGP is , where and , , with . Idiosyncratic terms are independently generated as with , for and defining the degree of serial correlation among factors, and with defined as diagonal and Toeplitz matrices. Common factors are estimated by principal components in cases with , and by the procedure of pena2006nonstationary in the case with . is the of a regression of actual on estimates factors. is the median of the MSE between the actual and estimated factor loadings. Haar and D8 wavelet functions are used in the study. All experiments are based on 1000 replications.

Table 1: , , and . The measure of the consistency of the estimated unobservable factors and the estimated factor loadings are presented in the report.

Table 1 shows the results of the estimations of (15) and (16). As can be seen, the methodology proposed in this paper performs very well in relatively small samples regardless of size distortion between and . As seen in Table 1, is relatively high, indicating the satisfactory performance of the estimator. However, precision is a bit reduced when increasing the value of . These findings are maintained for both types of wavelets used and even when we allow for cross-correlation between idiosyncratic errors. Furthermore, inspecting the MSE in Table 1, we find that the MSEs decrease as increases in all cases, even if the common factors are serially correlated or if idiosyncratic errors are cross-correlated. Furthermore, another finding indicates that, in general, factor loadings using the wavelet D8 perform better than the wavelet Haar. We think that such results are reasonable due to the smoothness of the wavelet D8 in contrast with the other one. The wavelet Haar should be implemented when factor loadings’ dynamic have breaks or perhaps some aggressive jumps. We do not go any further in this line. Wavelet Haar’s features are part of an another research and are out of the present scope.

Furthermore, in Figures 1 and 2, we display the methodology’s performance to estimate the factor loadings. We choose a couple of different smooth functions for each value of and in both types of wavelets for comparison purposes. We consider the following functions: i) , and ii) , where indicates the loading of the first factor of the cross-sectional unit , and , the loadings of the second factor of unit

. Figures display the actual and estimated time-varying loadings and their bootstrap confidence interval at 95% with

replications, following moura2012. In such figures, we can see that the methodology works well independently of the value taken in .

Finally, as a complement of this simulation study, we compare our methodology with the approach of mikkelsen2018consistent. Their estimation methodology is similar to ours in the first stage. The difference is in the second stage because instead of using wavelet functions, as in this paper, they employ Kalman filter procedures to estimate the likelihood function. It is worth mentioning that the main difference between both setups lies in the performance of loading factors. While we assume smooth variations in this paper, their approach consists of stationary VAR processes.

We use the same DGP as before, but we focus only on the following simplest case.

where the matrix is a diagonal matrix, and idiosyncratic terms are independent of each other. As before, the vector of loadings is generated alternately by some smooth functions. We split them into two groups; i) smooth sine/cosine functions and ii) smooth trending functions (linear, square root, exponential, and logs trends).

Our analysis finds that the methodology proposed by mikkelsen2018consistent does not adjust either smooth sinusoidal functions or smooth trending functions in most cases. The worst performance appears when loadings are trending functions, while it is possible to rescue from time to time some proper estimations when loadings behave as sine/cosine functions. Given our results, we argue that the model proposed in this paper seems to be more suitable when loadings have deterministic trends; in contrast, the methodology proposed by mikkelsen2018consistent fails to capture this type of behavior in loadings by treating them as a stationary VAR.

We repeated the entire simulation study as in Table 1 using the methodology of mikkelsen2018consistent; however, given the inferior performance of loading estimates, we prefer not to augment Table 1 for brevity. Instead, we support this part of the simulation study with Figure 1 that compares both approaches with four examples of smooth functions, two cases when loadings have trends, and two with sinusoidal functions.

Figure 1: Comparison between the methodology proposed in this paper (solid red line) and that by mikkelsen2018consistent (solid black line) to estimate time-varying loadings (dashed blue line) defined by the smooth function written above of each figure.

5 Application

This section applies our methodology to study comovements in loads and prices of the Nord Pool power market.

Nord Pool runs the leading power market in Europe and operates in the day-ahead and intraday markets. Elspot is the day-ahead auction market, where participants act in a double auction and submit their supply and demand orders (spot prices and loads) for each hour of the next day. The market is in equilibrium when demand and supply curves intersect at the system price and system load for each hour. The hourly system prices and loads series are announced as 24-dimensional vectors, which are determined simultaneously.

Electricity markets have particular features that do not exist in any other commodity market. Remarkably, the non-storability of electricity provokes that the time series of prices shows excessive volatility, possible negative prices, and many spikes over time. Other relevant features of electricity prices and loads are the intra-day, week, and year seasonal components, see weron2007modeling.

Univariate time series methods have mostly studied hourly electricity prices and loads, see Weron20141030 for a rich review. Some authors have explored these series by multivariate techniques as high dimensional factor models until the last years. Using data from the Iberian Electricity Market (MIBEL), seasonal factors have been extracted in the works of alonso2011seasonal and garcia2012forecasting, whereas alonso2016electricity propose to employ model averaging factor models to improve forecasting performance. Pennsylvania - New Jersey - Maryland (PJM) interconnection market is studied by maciejowska2015forecasting who estimate factor models for forecasting evaluation using hourly and zonal prices. Furthermore, the Nord Pool power market is studied in the works of ergemen2016common, who study the long-term relationship between system prices and loads, and rodriguez2017estimation, who use regional prices in a multi-level setting.

All these studies assume that factor loadings do not vary over time. In this sense, we are the first to study the power market with a factor model with time-varying loadings to the best of our knowledge. Such a setup may help to attract some dynamics not necessarily captured in common factors. We argue that factor loadings should evolve smoothly, perhaps reacting to the auction market’s smooth changes along the day and along year due to the series’ inherent seasonality.

We consider a balanced panel data set consisting of hourly prices and loads for each day from 13th March 2016 to 31st December 2018, yielding a total of daily observations in each hour. We download series from the Nord Pool ftp server, and prices are in Euros per Mwh of a load. Figures 3 and 3 display six time series in logs from which we can observe some characteristics in both time series. First, electricity system prices and loads vary differently along months with a common pattern in their evolutions. Second, the price series show many spikes which are related to the way how the market operates. Third, seasonal variations are much more marked in loads than in price series. Fourth, electricity prices and loads seem to have nonstationary behaviors, then, we set our estimation on the nonstationary approach discussed in section 2.2.

Figure 2: Hourly system loads in logs for six different hours showing working and non-working hours performances, 12 March 2016 to 31 December 2018.
Figure 3: Hourly system prices in logs for six different hours showing working and non-working hours performances, 12 March 2016 to 31 December 2018.
Figure 2: Hourly system loads in logs for six different hours showing working and non-working hours performances, 12 March 2016 to 31 December 2018.

We use the procedure proposed by alessi2010improved with first-differenced variables to estimate the number of factors. This methodology introduces a multiplicative tuning constant in the penalty function to improve the criteria of bai2002determining. In line with the literature, we find two common factors in prices, and the same number in loads. These factors explain 95% of the variation in the panel of prices, and 97% of loads. These percentages are slightly higher than those found in ergemen2016common, but our period is shorter and does not cover periods of infrastructure delineation in the Nord Pool power market. Therefore, we estimate the model in (12) with two common factors for the panel of loads and the panel of prices.

Figures 4 and 5 display estimates of common factors of hourly system loads and prices, respectively. As seen in Figure 4, the first factor of hourly loads captures the strong seasonal component showing possible weekly and monthly periodicities. The second factor seems to capture mainly a more erratic weekly variability, mostly during working hours, as seen in Figure 3. On the other hand, Figure 5 shows that estimated common factors of hourly system prices exhibit some stylized facts of the underlying series. Volatility clustering and nonstationary behavior are captured by the first factor, while excessive price spikes in 2018 by March, July, and August are extracted mostly by the second factor. ergemen2016common also document all these characteristics of common factors.

Figure 4: Common factors of hourly system loads.
Figure 5: Common factors of hourly system prices.

Figures 6 and 7 show time-varying loadings of the first and second common factors of hourly system loads, respectively. Figures 8 and 9 show the time-varying loadings estimates of hourly system prices. We only display the results for six hours for exposition purposes: 04:00, 08:00, 12:00, 16:00, 20:00, and 24:00 hrs. These hours help to observe different performances of loadings across working and non-working hours. For comparison purposes, we estimate the model proposed by mikkelsen2018consistent using the generalized covariance matrix introduced in (5), in the first stage, to allow for nonstationary variables. We standardized loadings estimates to compare both approaches.

Figure 6: Time-varying loadings of the first factor of hourly system loads. Black solid lines represents loadings estimated by the method proposed in this paper, while blue dot lines represents loadings estimated by mikkelsen2018consistent.
Figure 7: Time-varying loadings of the second factor of hourly system loads. Black solid lines represent loadings estimated by the method proposed in this paper, while blue dot lines represent loadings estimated by mikkelsen2018consistent.
Figure 8: Time-varying loadings of the first factor of hourly system prices. Black solid lines represent loadings estimated by the method proposed in this paper, while blue dot lines represent loadings estimated by mikkelsen2018consistent.
Figure 9: Time-varying loadings of the second factor of hourly system prices. Black solid lines represent loadings estimated by the method proposed in this paper, while blue dot lines represent loadings estimated by mikkelsen2018consistent.

Further inspection on Figure 6 indicates that factor loadings corresponding to the first common factor of system loads show smooth variations along time when estimating them by the model proposed in this paper (solid black lines in figures). These time-varying loadings capture a kind of readjustment of system loads when the demand reaches a maximum level (in winter) and a minimum level (in summer), which aligns with the market’s nature. Moreover, time-varying loadings in the remaining months oscillate smoothly around the mean. In contrast, factor loadings by mikkelsen2018consistent (blue dot lines) behave smoothly during non-working hours (see, 04:00, 20:00, and 24:00 hrs.), but the performance is very erratic during working hours (see, 08:00, 12:00, and 16:00 hrs.). Our intuition says that such erratic behavior comes from the volatile behavior of the hourly loads in the same working hours (see again Figure 3). Therefore, the main difference between both approaches appears in hours when series behave more volatile. However, it seems intuitive to think that factor loadings should behave smoothly throughout the year without considering the time of day because of changes in the electricity demand.

In turn, factor loadings corresponding to the second common factor of system loads also have a regular periodic behavior along time, see solid black lines Figure 7. It is interesting to see that loadings behave smoothly independently if the second common factor is much more volatile than the first one (see Figure 4). In contrast, the second approach (blue dot lines in the figure) reacts aggressively along time and could be related to the performance of the second common factor.

Concerning loadings estimates of hourly system prices using the methodology proposed in this paper (solid black lines in Figure 8), at first glance, unlike system loads, we do not find that factor loadings have precise periodic movements along time. In general, they also behave smoothly as before. It is worth mentioning that some moments where loadings seem to be time-invariant, see, for instance, the last year in Figure 6. Note that this behavior is not captured by loadings estimates using the approach of mikkelsen2018consistent (see blue dot lines in the same figures), which provides loadings with a very volatile performance.

This empirical study leaves many open possibilities for future research. First, testing if loadings are time-varying or constant over time. This test would improve the specification of factor models and help to understand the market behavior. Second, short-term forecasting of prices and loads is essential for energy markets. Therefore, another research may investigate if loads/prices forecasts obtained by time-varying factor models are more accurate than those provided by standard setups. Third, as discussed in the Monte Carlo study, the factor loadings using the wavelet D8 perform better than Haar in situations where time-series evolve smoothly, as in the case of system loads. In this sense, a deeper study on electricity prices could help us understand if the wavelet Haar represents a better choice for modeling the performance of factor loadings since electricity prices are more volatile than loads.

To conclude, our setup and methodology can be used in panel data models where unobservable common factors drive the cross-sectional dependence. Recently, rodriguez2020air explored this possibility to analyze the dynamic behind air pollution and mobility to disentangle the contagion in the COVID pandemic.

6 Concluding remarks

Factor models’ standard approaches assume that factor loadings are invariant along time. Here, we relaxed such an assumption allowing for time-varying loadings that behave as smooth and continuous-time functions. The paper is novel because time-varying can capture a more realistic behavior compared to some proposals in the literature that assume that stable stochastic processes drive loadings. Our estimation methodology is a two-step procedure based on GLS with wavelet functions. We explored stationary and nonstationary setups to the extent of the model’s applicability. We recommend using Wavelet D8 estimates in empirical applications with smoothed loadings based on our finite sample analysis. In this sense, Wavelet D8 provides better estimations than Wavelet Haar when factor loadings do not have sudden changes.

Finally, in our empirical study, we focus on the complex dynamics of Nord Pool electricity loads and prices in a large panel of hourly observations. Our findings indicate that factor loadings vary smoothly over time in prices and loads. Notably, in loads, the time-varying loadings have a periodic behavior, which obeys the market’s seasonality. The analysis provides relevant insights into the dynamics of the market. Future research will point to how the model proposed can be used from a forecasting perspective.


Appendix A Technical appendix

a.1 Proof of Theorem 1