share
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
| (1) | |||||
| (2) |
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:
-
For , there is a function with
such that
with and is a positive definite diagonal matrix.
-
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.
Factors:
-
[label=]
-
-
as with is a positive definite diagonal matrix.
Assumption B.
Factor loadings:
-
[label=]
-
and , as , with positive definite matrix , and is the th row of
-
where is the th row of
-
, for and
Assumption C.
Idiosyncratic terms:
-
[label=]
-
-
with , for a constant and
-
and
Assumption D.
Time-varying factor loadings and factors:
, and are functions such that the following conditions are fulfilled for any and
-
[label=]
-
,
-
-
.
Assumption E.
Independence:
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
andensure 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
| (3) |
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
| (4) |
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
| (5) |
where and can be either 0 or 1.
Under this framework, consistent factor estimates are given by
| (6) |
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
| (7) |
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
| (8) |
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
| (9) |
with , and In matrix form, we have
| (10) |
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 ,
| (11) |
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
| (12) |
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
| (13) |
where
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 bywhere 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
where
- 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
| (14) |
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:
-
The precision of the estimation factors is measured by the statistics as in bates2013consistent, given by
(15) 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.
-
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
(16)
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 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.
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.
 |
 |
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.
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.
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.
References
Appendix A Technical appendix
a.1 Proof of Theorem 1
Proof.
Let