Periodic dynamic factor models: estimation approaches and applications ∗

: A periodic dynamic factor model (PDFM) is introduced as a dynamic factor modeling approach to multivariate time series data exhibiting cyclical behavior and, in particular, periodic dependence structure. In the PDFM, the loading matrices are allowed to depend on the “season” and the factors are assumed to follow a periodic vector autoregressive (PVAR) model. Estimation of the loading matrices and the underlying PVAR model is studied. A simulation study is presented to assess the performance of the introduced estimation procedures, and applications to several real data sets are provided.


Introduction
Dynamic factor modeling of multivariate (vector) time series is a common and popular approach, especially when the dimension of a time series is large, possibly larger than the sample size. See, for example, Bai and Ng (2008), Stock and Watson (2011), Reichlin (2011, 2012). A dynamic factor model (DFM, for short) of a q-vector time series {X n } n∈Z is expressed as where μ is the mean of X n , Λ 0 is a loading matrix of dimension q × r, {F n } is a r-vector time series consisting of r factors with r typically much smaller than q and {ε n } is a q-vector time series uncorrelated with {F n }, sometimes assumed to be white noise (i.e. {ε n } ∼ WN(0, Σ ε )) though it can also be correlated in time. The factors {F n } are commonly assumed to follow a vector autoregressive model of order p (VAR(p) model), where Φ(B) = I − Φ 1 B − . . . − Φ p B p is a VAR matrix polynomial with the usual backward shift operator B, each matrix Φ i of dimension r × r, and {ζ n } ∼ WN(0, Σ ζ ). In fact, the DFM (1.1)-(1.2) is also called static; more generally, DFM is defined by replacing Λ 0 in (1.1) by a loading matrix polynomial Λ(B). Depending on the assumptions and the estimation method used, additional restrictions on the model parameters are also made (e.g. Λ 0 is lower triangular, Σ ε is diagonal and Σ ζ = I r , the r × r identity matrix). The model (1.2) is also commonly assumed to be strictly stationary and causal, though nonstationary DFMs have also been considered.
In this work, we are interested broadly in dynamic factor models for multivariate series exhibiting cyclical behavior with period s. We use the term "cyclical behavior" loosely and rather mean by it a collection of associated models. Two main modeling approaches for univariate series exhibiting cyclical behavior are based on either (multiplicative) seasonal ARIMA models (e.g. Brockwell and Davis (2009), Ghysels and Osborn (2001)) or periodic ARMA models (e.g. Franses and Paap (2004), Hurd and Miamee (2007)). Dynamic factor modeling based on seasonal vector ARIMA models has been studied by several authors. Alonso, García-Martos, Rodríguez and Sánchez (2011) assume that common factors follow a seasonal VARIMA(p, d, q)×(P, D, Q) s model where s is the period, Φ(B) is as in (1.2), Φ s (B s ) = I − Φ s,1 B s − . . . − Φ s,P B sP is a seasonal VAR matrix polynomial, Θ(B) and Θ s (B s ) are, similarly, MA matrix polynomials of respective orders q and Q, 1 and d, D are non-negative integers (that allow for factors to be non-stationary when d ≥ 1 or D ≥ 1). Using a state space formulation of the model (1.1)-(1.2), and imposing flexible restrictions on model parameters, the authors developed an EM algorithm for model estimation. Furthermore, bootstrap strategies were proposed to quantify statistical uncertainty associated with parameter estimates and forecasts. Nieto, Pena and Saboyá (2016) proposed a similar model but identifying separately the factors that are (i) non-stationary and non-seasonal, (ii) (non-stationary) seasonal, and (iii) stationary. A notable procedure to estimate the dimension of the three different types of factors was proposed by the authors, based solely on suitable covariance matrices. In contrast to seasonal models, there has seemingly been no work on dynamic factor modeling of multivariate series for which periodic models are better or more natural alternatives than the seasonal models above. The goal of this paper is to introduce and to study several approaches to such periodic dynamic factor modeling. The most general form might be to assume, in the "stationary" and "static" form, that X n − μ m = Λ m F n + ε n , (1.4) where m = m(n) = 1, . . . , s, refers to the season that n belongs to, the q × r loading factors Λ m possibly depend on the season m, {ε n } is a WN series whose covariance matrix Eε n ε n = Σ ε,m possibly depends on the season m (denoted {ε n } ∼ WN(0, Σ ε,m )) and the r ×1 factors {F n }, which are independent of {ε n }, are assumed to follow a periodic VAR model of order p (PVAR(p), for short), (The order p may depend on the season as well but we consider the case of a single p for simplicity.) We shall refer to (1.4)-(1.5) as a periodic dynamic factor model (PDFM). As with (1.2), we also assume that (1.5) is in the stationary and causal regime. We note that as defined above, the PDFM is quite general. It includes several potentially interesting special cases including: • periodic Λ m , Σ ε,m , non-periodic stationary F n ; • constant Λ m = Λ, Σ ε,m = Σ ε , periodic F n ; (1.6) • periodic Σ ε,m , constant Λ m = Λ, non-periodic stationary F n ; and others (but see also Remark 2.1 below). In all these special and general cases, the component series of the PDFM (1.4)-(1.5) are periodically correlated in the sense that the mean μ(n) = EX n and the autocovariance function γ(n 1 , n 2 ) = Cov(X n1 , X n2 ) are periodic with period s, that is, μ(n) = μ(n + s), γ(n 1 , n 2 ) = γ(n 1 + s, n 2 + s), n,n 1 , n 2 ∈ Z. (1.7) Hence, PDFMs should be used primarily when periodic correlations are observed for univariate component series. This will be illustrated on real time series below. In this work, we shall study estimation methods for the PDFM (1.4)-(1.5) and some of its special cases as in (1.6). Following the approaches taken for DFMs in the literature, we shall consider two estimation frameworks for the so-called exact and approximate model formulations. (The distinction is made and discussed nicely in Doz et al. (2012).) For the exact DFM formulation, it is assumed for identification purposes that Σ ζ is the identity matrix and Λ 0 is lower triangular, and without sacrificing much model flexibility, it is further assumed that Σ ε is diagonal (see, e.g. Peña and Poncela (2006), Alonso et al. (2011)). Sometimes, the polynomials Φ(B) are assumed to be diagonal as well (e.g. Nieto et al. (2016)). Other identifiable model formulations are also used (e.g. Bai and Wang (2015)). For the exact model formulation, estimation is typically carried out using the EM algorithm. We shall also make similar assumptions on our PDFM in the exact case and employ the EM algorithm. The exact model though is usually fitted for smaller dimensions q, which we shall also do here.
In higher dimension q, assuming that both q and T , the sample size, converge to infinity, it has been well recognized that weaker assumptions on the model factor structure can be imposed, under the term "approximate DFM", which lead to consistency results subject to properly rotated factors. See Doz et al. (2011Doz et al. ( , 2012, Ng (2007), Forni, Hallin, Lippi andReichlin (2000). Moreover, in this case, the model parameters in approximate DFMs can be estimated using principal components, though the effects of additional iteration steps involving smoothing as in the EM algorithm have also been studied. In this work, we shall extend the approximate DFM formulation and convergence results to PDFMs, relying heavily on the work by Doz et al. (2011Doz et al. ( , 2012). An additional important feature of this work is that we shall also be using sparse models on multiple occasions. This should come as no surprise since allowing DFM coefficients to depend on the season increases their number considerably.
The rest of the paper is organized as follows. Estimation of PDFMs is considered in Section 2, including the statements of the main results of this work. A competing approach where standard DFM is applied to the multivariate series obtained by "stacking" X n 's across the s seasons, is also discussed in Section 2. A simulation study is presented in Section 3, and applications to several real data sets can be found in Section 4. Section 5 concludes while Appendix A includes the proofs of the results mentioned in Section 2, and Appendix B discusses the adaptive lasso estimation procedure for the so-called Fourier PVAR model.

Estimation of PDFMs
We present here estimation methods for PDFMs in the two formulations, exact (Section 2.1) and approximate (Section 2.2), as discussed in the context of DFMs in Section 1. Some asymptotic results for the approximate formulation are also given (Section 2.3) along with estimation methods based on sparsity and the Fourier representation (Section 2.4). The case of a different number of factors across seasons in (1.4) is also discussed (Section 2.5).

Exact model formulation and EM estimation
In the exact formation, we assume that the PDFM parameters in (1.4)-(1.5) satisfy: where B m has dimension (q − r) × r. For the usual DFMs, the identifiability of model parameters under the assumption analogous to (2.1) was proved in Bai and Wang (2015), and their arguments extend straightforwardly to PDFMs under (2.1). In fact, a basic reason for (2.1) leading to identifiability is that essentially any q × r loading matrix Λ m can be reduced to the form (2.1) as where Λ (1) m and Λ (2) m are r × r and (q − r) × r, respectively, and the r × r matrix C m is such that Λ (1) m C m = I r (which requires that Λ (1) m is non-singular). Since Λ m F n = Λ m C m C −1 m F n , the use of C m to achieve the identifiable form (2.2) just means that the factors F n similarly need to be adjusted through C −1 m . As in the usual exact DFMs, the parameters of the models (1.4)-(1.5) can be estimated with the EM algorithm. For this purpose and e.g. when p = 1, we have adapted here the arguments in the R code for the EM algorithm found in Shumway and Stoffer (2006), pp. 342-344, with the following changes: in the M-step of the algorithm, we also need to minimize with respect to the loading matrix (A t in Shumway and Stoffer (2006)); we allow both the loading matrix A t and the PVAR matrix Φ t to depend on the season, and, finally, in the M-step, we also work and perform minimization under the parametrization (2.1) of the loadings. The R code implementing this will be available online upon publication of the article. We also note that the EM algorithm requires initial values of the model parameters. For these, we use the values obtained from the principal component approach discussed in Section 2.2 below, after transforming the latter to the used identifiable parametrization through (2.2). As with all likelihood calculations, the EM algorithm can also yield confidence intervals for the parameters of interest, through inverting the numeric Hessian.

Approximate model formulation and PCA estimation
We discuss here estimation of factors and the underlying PVAR model in the PDFM defined by (1.4)-(1.5) in its approximate formulation. As noted in Section 1, the latter refers to model assumptions and the fact that large q is considered.
One of the basic observations is that PDFM can be viewed as consisting of s separate DFMs. That is, let X ( 2.3) The factors {F (m) t } can then be estimated as in the usual DFM. The exact procedure is adapted from Doz et al. (2011), and is recalled below. Estimation of the underlying PVAR model in (1.5), which is also discussed below, requires a more delicate extension. As in Doz et al. (2011), we assume without loss of generality that EF and Λ m for season m = 1, . . . , s, are given by where r is the number of factors used.
4. Use the least squares method to estimate the parameters of the PVAR(p) model in (1.5) based on F n , n = 1, . . . , T s, to obtain the estimates Φ m,1 , . . . , Φ m,p and Σ ζ,m , m = 1, . . . , s. The covariance matrix Σ ε,m = Eε Several comments regarding the steps above are in order. The estimators (2.6) are defined as they appear in the literature of principal component analysis, e.g., Bai and Ng (2008), Stock and Watson (2011). The principal component estimators in (2.6) of Step 3 are the solutions F (m) t and Λ m minimizing the sum of squares (2.7) subject to the normalization q −1 Λ m Λ m = I r . The last step, Step 5, follows Doz et al. (2011). We also note that iterating Steps 4 and 5 is possible and is closely related to employing the EM algorithm. That is, based on the final estimate F n of the factors, we can get new estimates of PVAR(p) parameters and minimizing (2.7) for given F n subject to normalization on loadings gives updated estimates of loadings, and so on. This is essentially the same as running the EM algorithm of Section 2.1, as used by Engle and Watson (1981) in the case of DFMs. Doz et al. (2011) also mention iteration procedures but their simulation study shows that just one iteration performs reasonably well in practice. Our simulation study comparing the exact and approximate model approaches in Section 3.2 also shows that just one iteration seems to work well in practice. The expression in (2.5) follows the convention used by Stock and Watson (2011). Finally, we should also note that the "decoupling" of the original PDFM (1.4)-(1.5) into s separate DFMs should not be too surprising. For example, estimation of time-varying loadings by means of a time-localizing kernel function can be found in Su and Wang (2017). As noted above, it is a bit more delicate to show that the estimated factors can be "pieced" back together into a PVAR model (1.5), assuming that the latter holds at the PDFM level.
Remark 2.1. Related to the last point, we also make the following important, even if seemingly trivial, observation. The factors F n in (1.4)-(1.5) can be identified up to rotations only, which further can depend on m. In this regard, note that the PVAR model (1.5) is "closed" under the multiplication of F n by such rotation matrices, in the sense that this would only change the parameters of the PVAR model. Indeed, with rotations C m and e.g. when p = 1, we have that Thus, the new rotated factor process C m(n) F n is necessarily a PVAR. This precludes the first case in (1.6) as a possibility in the model specification (1.4)-(1.5), and this case should not be considered.

Asymptotic results for approximate model formulation
The results proved in Doz et al. (2011) can also be used directly to describe the asymptotics of the estimated factors and the PVAR model parameters in the approximate model formulation. This shall require q, T → ∞; otherwise, no consistency can be expected. We shall suppose that the number of factors r and the underlying PVAR model order p are known -their estimation is discussed below. For identifiability purposes, we also need to fix a rotation up to which the factors and the PVAR model parameters will be identified. Thus, where D m is a diagonal matrix consisting of the eigenvalues of Λ m Λ m and Q m is the orthogonal matrix consisting of the corresponding eigenvectors. Set so that the PDFM can be expressed as ). Let also D m,r = diag( d m,1 , . . . , d m,r ) and P m,r = ( U m,1 , . . . , U m,r ) with d m,i and U m,i defined following (2.5) and set The estimates of Ψ m (B) based on G n will be denoted Ψ m (B). Note that F As in Step 5 above, we let G n be the smoothed estimate of G n obtained through the Kalman smoother.
The following assumptions will be made for the PDFM (1.4)-(1.5) following Doz et al. (2011). Basically, the assumptions (A1)-(A3) impose a factor structure characterized by stationarity and weak temporal dependence of the common and idiosyncratic components. The summability of the coefficient matrices in the Wold decompositions in (A3), in particular, excludes long-memory series. We note that the moment conditions across component indices i in (A2) and (A3) are assumed to hold irrespective of q → ∞. The condition (A3 ) posits a parametric PVAR(p) model of the common factors.
We shall further assume weak correlations on idiosyncratic components and relatively pronounced eigenvalues of the underlying factor structure in the sense of approximate factor model of Chamberlain and Rothschild (1983). Let A = (λ max (A A)) 1/2 be the spectral norm of a matrix A with λ max (A A) being the largest eigenvalue of a matrix A A. Let also λ min (A A) be the smallest eigenvalue. For any m = 1, . . . , s, assume:  Wang and Fan (2017). It remains to be seen whether their results can be extended to the dynamic context.
The first result concerns the asymptotics of the estimated quantities in the relation (2.10).
Remark 2.3. We expect that the result analogous to (2.14) holds for the smoothed estimated factors G n as in Doz et al. (2011). For brevity, we shall not pursue this result here -we found the effect of smoothing to be of little practical significance.
Several methods have been proposed to estimate the number of dynamic factors r. For a DFM, one popular approach is based on the information criteria studied by Bai and Ng (2007). We adapt this information criteria to find the number of factors for PDFM by setting ( 2.18) where SSE(k) is the sum of squared errors assuming k factors calculated by and g(q, N ) is one of the following four functions: where q ∧ N = min(q, N ). We also mention a popular graphical method of selecting the number of factors through the so-called scree plot. Adapting it to our context, the eigenvalues of the covariance matrix Σ Y from the deseasonalized series X n − μ m are plotted on the y-axis from the largest to the smallest, and the leveling-off point is selected as the number of factors r. See Figure 8 below for an example. A further discussion related to the choice of r can be found in Section 2.5 below. The order p of the underlying PVAR model is also estimated using an information criteria, namely, the BIC, where Σ ζ,m (k) is the covariance matrix estimate of PVAR(k) model residuals at season m.

Other approaches based on sparsity and DFMs
In the simulation study presented in Section 3 and the applications presented in Section 4, we also estimate the underlying PVAR model by using several other estimation methods, than just a simple least squares regression. On the one hand, we shall use sparse estimation, following our work in Baek, Davis and Pipiras (2017). Sparsity is natural and convenient to use with PVAR models, which have a large number of parameters to estimate even for moderate values of dimension r. On the other hand, we shall also use sparse estimation but in the Fourier representation of PVAR coefficients. Representing the univari-ate PAR coefficients in the Fourier basis is used in PAR model estimation in Jones and Brelsford (1967), Dudek, Hurd and Wójtowicz (2016a) (see also references therein), with the idea that only several of the Fourier coefficients will be needed assuming the PAR coefficients change slowly over a period. Extending the Fourier representation approach to PVAR models, the (j, k) element of Φ m,i is represented as (2.22) where H is the order of the Fourier approximation. For example, if r = 5, s = 24 and p = 2, then the number of parameters for PVAR(2) model is given by r 2 ps = 1200. While in the Fourier approximation with H = 2, the number of parameters to be estimated is reduced to r 2 p(2H + 1) = 250. We can further impose sparsity on the Fourier coefficients to estimate them using adaptive lasso as described in Appendix B. We shall refer to the modeling approach based on (2.22) and adaptive lasso as Fourier-PVAR or FPVAR, for short.
Further insight into PDFM can be given through the following discussion. PVAR time series {F n } satisfying (1.5) are nonstationary in general but the multivariate time series F t = vec(F (1) t , . . . , F (s) t ) obtained by stacking F n 's over one full period is stationary (see also the notation around (2.3)). In fact, the series { F t } has a VAR representation, and this is why additional suitable assumptions on the stationary "block" series {F (m) t } can be made (as in Section 2.3). These observations also naturally suggest the following "competing" dynamic factor approach to multivariate time series with periodic dependence structure. As with the PVAR model above, a natural possibility is to consider the "block" vector series X t = vec(X where the various components of the model are interpreted as in (1.1)-(1.2). This "stacking" approach was used by e.g., Dordonnat, Koopman and Ooms (2012), in the case of univariate series X n . For later use, we shall refer to the approach based on (2.23)-(2.24) as block-DFM (where "block" refers to the blocks X (m) t used to construct X t ). Note that block-DFM becomes our PDFM when Λ 0 and Φ(B) have block-diagonal structures, namely, (2.25) That is, PDFM can be viewed as block-DFM with restrictions (2.25). This means, in particular, that a larger number of parameters needs to be estimated in block-DFM: the loading matrix Λ 0 has dimension sq × sr and hence s 2 qr parameters, and the polynomial Φ(B) has p(sr) 2 parameters, as opposed to the respective numbers sqr and spr 2 for PDFM. As a consequence, if the data actually follows PDFM, we can naturally expect better estimation performance when going the PDFM rather than the block-DFM route. A smaller number of parameters should also lead to better forecasts -incidentally, this is also what we generally observe in data applications (see Section 4 below).
The number of parameters aside, the option of using PDFMs versus block-DFMs is similar to that of using PVAR models versus VARs (or larger dimensional embedded VARs). Discussions for the latter case can be found in Hurd and Miamee (2007) p. 6, Franses and Paap (2004), Section 1.3. Among advantages of periodic models, their flexibility and convenient interpretation are often mentioned. In some applications, periodic dependence structures are also naturally expected: e.g., in environmental applications, this structure may be determined by the annual seasons, or in engineering applications, as in vibrations of heavy machinery, they may reflect the periods associated with machinery rotating parts.

Varying number of factors
We assumed in the preceding sections that the number of factors is the same across the s seasons. But the approaches described in Sections 2.1 and 2.2 extend straightforwardly to the case of varying number of factors, but with several caveats that are mentioned below. We do not consider this more general case throughout the paper for two reasons: simplicity and also since the real data considered in this work do not present strong evidence in favor of using varying number of factors.
To be more specific, by the case of varying number of factors, we mean the model (1.4)-(1.5) but with where Λ 0,m has dimension q × r m , 0 q×(r−rm) is the q × (r − r m ) matrix of zeros, r = max m r m and at least one of r m is smaller than r. In other words, only the first r m components of the r-vector factors F n load onto X n . Consider first such a model in the framework of Sections 2.2-2.3, where q is assumed to be large. Since Steps 1-3 of the estimation algorithm in Section 2.2 apply for fixed season m, they can also be used for the model (1.4)-(1.5) with (2.26). Moreover, a different number of factors might be suggested by using the criterion (2.18) but replacing SSE(k) in (2.19) by SSE m (k) in (2.7), or also by examining the scree plots of the eigenvalues of the matrices Σ (m) Y . The factors that load in season m can still be estimated consistently by the results analogous to those in Section 2.3.
Step 4 of that algorithm, however, does not apply. One is presented with an interesting problem of estimating the parameters of a PVAR model but where for some season(s) m, only r m -vector parts of the r-vector series F n are actually observed. This issue of "missing" data could be dealt with naturally through the EM algorithm as e.g. in Shumway and Stoffer (2006), Section 6.4. But one question that remains to be addressed (and that we currently do not know how to answer) concerns identifiability issues for such PVAR models with "missing" data. What we mean can be explained through the following simple example. Suppose F n is one-dimensional and satisfies a PAR(1) model with period s = 4 and the AR parameters φ 1 , φ 2 , φ 3 and φ 4 for the four seasons. If, for example, F n is "missing" in the periods 2 and 3, it is unrealistic that φ 2 and φ 3 can be identified since there is no data to deduce the chronological order of seasons 2 and 3 (that is, the model parameters φ 2 and φ 3 may as well be interchanged).
A similar discussion also applies to the model (1.4)-(1.5) with (2.26), in the framework of Section 2.1 where q is thought to be small. That is, modulo the identifiability question raised above, there are a priori no major issues of running the EM algorithm for the model (1.4)-(1.5) with (2.26) in its state space formulation. One smaller issue to mention is that in the M-step of the algorithm, we do not have a closed-form solution to Λ m and use direct optimization. The EM algorithm was tested successfully on simple models with varying number of factors, and the associated code is available upon request.

Simulation study
In this section, we study the finite sample performance of the proposed estimation methods for PDFMs by Monte Carlo simulations. We focus on the PCA method of Sections 2.2-2.3 and then present briefly the results for the EM algorithm of Section 2.1.

Approximate model formulation with PCA
For brevity, we consider loadings and idiosyncratic error terms that do not depend on the season (but naturally the factors follow a PVAR model). The univariate components X j,n of X n in the PDFM (1.4) can then be written as X jn = r k=1 λ jk F kn + ε jn , j = 1, . . . , q.
The factor loadings λ jk of the matrix Λ = (λ j,k ) are generated as the (j, k) components of the QR decomposition of q × r matrix consisting of i.i.d. N (0, 1) random variables and multiplied by √ q so that 1 q Λ Λ = I r (see, in particular, Assumptions (CR1) and (CR2) in Appendix A). The idiosyncratic component is assumed to follow the AR(1) model ε jn = ρε j,n−1 + ξ jn (3.1) with ξ jn being i.i.d. N (0, 1 − ρ 2 ) so that this yields serial correlations with zero mean and unit variance. Regarding factors F kn , we considered two PVAR(1) models. The first model (DGP1) parameters are given by: for s = 4 and r = 3, The second model (DGP2) parameters are: for s = 2 and r = 4, Since factors are identifiable up to a rotation (see Section 2 and also Remark A.1 in Appendix A), we need to reestimate them and PVAR coefficients by fixing a particular rotation. We estimate the rotation matrix Q in (2.8) (which does not depend on m by our assumptions) as in Doz et al. (2011), p. 195, Then, the estimates of the PDFM factors and loadings are defined asF n = Q G n , Λ = Π m Q , and the PVAR model parameters Φ m,i are re-estimated based on F n , denotedΦ m,1 , m = 1, . . . , s.
We measure the performance of the estimation procedure through the following statistics. The first statistic measures the discrepancy of the estimated factors through a trace R 2 of the multivariate regression ofF onto F . Introduced in Stock and Watson (2002), it is defined as Table 1 The performance of the PDFM for DGP1 measured by R 2 and χ 2 (F , F ).  Table 3 The performance of the PDFM for DGP2 measured by R 2 and χ 2 (F , F ). where E denotes the empirical expectation calculated by taking average over Monte Carlo replications. The second summary statistic is a chi-square type statistic for common component used in Forni et al. (2000). It is given by The last summary statistic is a mean squared error of the loading matrices and PVAR parameters given by We considered three different dimensions q = 20, 50 and 100 with sample sizes T = 100 and T = 200. For the idiosyncratic errors following AR(1) model, we considered ρ = .3 and .7. All results are based on 500 Monte Carlo replications. Tables 1 and 2 report on the performance of PDFM for DGP1. It is observed that R 2 approaches 1 and χ 2 (F , F ) diminishes as the sample size T Table 4 The MSEs of the PDFM model for DGP2. increases and the dimension q increases, as in Doz et al. (2011). This means that the precision of the factor estimation increases as the data dimension increases in both directions. However, the MSE of loading matrices are increasing as the dimension q increases while the sample size T is fixed. The MSEs for PVAR coefficients are all decreasing as the sample size T increases, and sparse PVAR (sPVAR) method shows the smallest MSEs as we expected since the true model is sparse. Note also that the performance of PDFM is getting worse under stronger correlations for idiosyncratic errors. Similar results are observed for DGP2 as shown in Tables 3-4. It is interesting that Fourier-PVAR (FPVAR) model gives the smallest MSEs amongst all the methods considered in this particular simulation setting. This is because, as can be checked, the number of zero coefficients in the Fourier representation is even larger than that in the Φ m -representation of the DGP2 model. The average estimated PVAR coefficients are illustrated in Figures 1-2 with the coloring and its lightness associated with the displayed numerical values. In summary, the proposed methods for estimating PDFM perform well in finite sample.

Exact model formulation with EM
We report here briefly on the performance of the estimation procedure in the exact model formulation with the EM algorithm introduced in Section 2.1. The data generating process (DGP3) uses the identifiable model assuming (2.1). We set q = 4 with r = 2 and s = 4 with common loadings over the seasons. The underlying factors follows the PVAR(1) model with and idiosyncratic component follows the AR(1) model (3.1) with ρ = .3 or ρ = .7. We considered T = 100 and T = 200 and the performance measures are MSE(Λ, Λ) and MSE(Φ, Φ). The initial estimates for the EM algorithm are based on PCA, that is, the principal component estimator (2.6) is rotated as in the relation (2.2).   Table 5 summarizes the results. We compared the EM results with the twostage PCA estimator. For fair comparison, we rotated the estimates from the two-stage estimation by multiplying C based on the relation (2.2) and reestimated the PVAR parameters from the rotated factors C −1 F n . Observe first that the MSE decreases as the sample size T increases and the autocorrelation in idiosyncratic component is getting weaker. In short, the EM algorithm successfully estimates the underlying parameters. But note also that our two-stage PCA estimator is comparable to EM even if the dimension q of the data is 4, that is, quite low.

Applications
We apply here PDFMs to several real data sets. As noted in Section 1, PDFM should be entertained only in the case when univariate component series exhibit periodic correlations. In the applications below, the latter is the first thing examined and discussed. We also note that the detection tools not only reveal periodic correlations but also suggest what the period s is. When fitting a PDFM, the number of factors needs to be estimated as discussed at the end of Section 2.3. This is also one of the first tasks carried out in the applications below. The rest of PDFM analysis involves estimation of the loadings Λ m and the parametric model of the factors F n . Here, we shall consider a number of possibilities, discussed in more detail below, with these model components being either periodic or non-periodic, with estimation being either sparse or not, etc. Furthermore, the results will be compared to a number of alternative approaches, for example, the most elementary one where just the univariate series are being modeled.

Fraser river water flow
The first data set considered here is a monthly water flow from the Fraser River in British Columbia (BC), Canada. The Fraser River is the longest river in BC and the largest river by volume flowing into the Pacific seaside. It runs almost 1400 km and the annual discharge at its mouth is 112 cubic kilometers. The analyzed data are the monthly averages of the daily discharge measurements, in cubic kilometers. Six locations are selected: the Fraser River above Texas Creek, at Hansard, at McBride, at Hope, at Red Pass and at Shelly from 1959 to 2010; hence 624 observations are used in the analysis. We also reserved the last 24 observations for out-of-sample forecasts. The data set is obtained from the Wateroffice of Canada -see http://wateroffice.ec.gc.ca for more information.  The data is log-transformed and Figure 3 shows time plots of log-transformed monthly water flow of Fraser river at the six locations. It is clearly observed that the data show cyclic variations. Figure 4 shows seasonal mean and standard deviation at the Fraser above Texas creek. The water flow seems low during the winter season, from December to May, then increasing rapidly from April, achieving the highest flow in July, and then slowly decreasing afterwards. The standard deviation is "bimodal" having peaks at April and October, so the water flow dynamics changes around those times. That is, the water flow dynamics enters a higher water flow period from May to September, and then changes to  a slow period from November to March. Figure 5 presents two commonly used plots to detect periodic correlations in univariate series (after removing its seasonal mean): the spectral coherence plot according to Hurd and Gerr (2008), and a related test statistic with a critical value line from Lund (2011). The spectral coherence is plotted using the R package perARMA (Dudek, Hurd and Wojtowicz (2016b)). The two plots are for the univariate series at Hope. If a series exhibits periodic correlations at period s (after removing the seasonal mean), the spectral coherence plot should have diagonal lines emerging at multiples of the index N/s, here N/s = 600/12 = 50. The plot in Figure 5 indeed suggests that the first major diagonal line is around the index 50. This could be determined easier from the Lund test statistic plot, which essentially averages the spectral coherence statistic at different indices along the corresponding diagonals, and also provides a critical values as the dashed line. Note that the Lund test statistic value exceed the critical value at the indices 50 and 100, suggesting that periodic correlations are present in this univariate series. Figure 6 presents the sample cross-correlations between some selected locations prewhitened by fitting univariate AR(1) models. All four panels show significant lag 1 (cross-) correlation and suggest that the multivariate modeling is also potentially necessary to account for between-location dependence.
To account for such periodically correlated observations, we first applied sP-VAR (sparse PVAR) models as in Baek et al. (2017). The best model is sP-VAR(1;110), where 110 indicates the number of non-zero coefficients, by minimizing the BIC and using adaptive lasso to estimate coefficients. Figure 7 shows the autoregressive coefficient matrix stacked by 12 months. For example, the first 6 × 6 matrix on the left panel shows the autoregressive coefficient matrix for January, the 6 × 6 submatrix below corresponds to February, and so on. It is observed that coefficients from December to March look similar by having larger coefficients on the diagonal and little-to-none correlation for the off-diagonal entries. The coefficients for April are similar to December-March, but having smaller coefficients and negative coefficients appear first at the Shelley location. In May, the diagonal entries shrink; then negative correlations amongst different locations are prominent in June and July. Then, cross-correlations are dramatically weakened in August and within series correlations become dominant with a few exceptions through September to November. This shows that the water  flow in the Fraser goes through several changes over a year. That is, during the winter period, December to March, the water flows appear almost decorrelated across locations but having strong within-series dependence. From April, within-series dependence is getting weaker, but between-locations dependence is getting stronger and peaks during June and July, and then rapidly disappears. This is consistent with the observations from Figure 4.
The results for PDFM are the following. Using the scree plot and the information criterion, the number of factors is selected as r = 2; so the PDFM reduces dimension from 6 to 2. Figure 8 shows the scree plot and time plots of the two factors found in the PDFM (two factors are also commonly selected when we considered scree plot for each season).
Twelve loading coefficients corresponding to each season are overlaid in Figure 9. Observe that there are some similarities in the shapes of loading coefficients, and if these are grouped by taking averages, Figure 10 is obtained. The grouping follows three patterns: one from December to March, then April to November except May, and May. Similar to the earlier interpretations, there is  a seasonal shift in water flow in the Fraser, mainly during the winter season and the rest of time. However, a further finding is that the water flow dynamics in May really differs from that for other months. It almost moves in the opposite direction and can be interpreted in that the negative between-locations correlations are dominant at that time.
To evaluate the forecasting ability, we also obtained the h-step-ahead out-ofsample forecasts and computed the empirical mean squared error (MSE), where X n+h is the h-step-ahead best linear predictor of X n+h based on the training sample size N − N r and the test sample size N r − h + 1. Table 6 shows the h-step-ahead out-of-sample forecast error together with the number of non-zero parameters for h = 1, 2 and 4. In the column "Models", DFM corresponds to a dynamic factor model with common loadings across seasons while DFMS stands for DFM with seasonal loadings. For example, DFM-PVAR is the model where the common loadings are assumed for the dynamic factors and PVAR model is applied for the estimated factors. The DFMG represents DFM with grouped loadings as in Figure 10. As before, sPVAR corresponds to sparse PVAR obtained using adaptive lasso, FPVAR to Fourier-PVAR, and PAR to univariate periodic autoregressive model applied to each component series. EM-DFM-PVAR refers to DFM-PVAR estimated using the EM algorithm as discussed in Section 2.1; all other dynamic factor models in the table use PCA estimation.
The results are mixed and no method outperforms in all cases. For the 1step-ahead forecasting, FPVAR(1) model shows the smallest MSE while DFMS-FPVAR(1) performs best for h = 2 and 4, and DFM-FPVAR(1) is also quite comparable, while using 6 fewer parameters than FPVAR(1). In general, it is observed that the PDFM with FPVAR show better forecasting performance than full models, while using fewer parameters. As seen from Figures 8-10 and the discussion above, the PDFMs also have nice interpretations. EM-PDFM-PVAR(1), which is based on the estimation method suggested in Section 2.1, does not assume sparsity and generally performs worse than its sparse analogues; but note also that it outperforms (in almost all cases) the non-sparse DFM-VAR(1), DFM-PVAR(1), DFMS-PVAR(1) and is comparable to the sparse DFM-sVAR(1) (their performance measures only differ in the fourth decimal places).
We also report here briefly on comparison to forecasting using block-DFM (see Section 2.4). Since the block-DFM takes the full period as one block, it is not fair to compare one-step-ahead out-of-sample forecasts. Instead, we calculated the sum of 12-step-ahead forecasting errors, where the sum of m-step-ahead forecasting errors is given by and X N −m (h) is the best linear predictor of X N −m+h from the training sample {X 1 , . . . , X N −m }. The block-DFM performs worse than the PDFMs, for example, the forecasting error is 153.37 versus 131.70 for DFM-sPVAR.

US quarterly macroeconomic data
We also considered the quarterly macroeconomic time series of e.g., GDP, GNP and industrial product index for the United States analyzed in Stock and Watson (2009 We checked for periodic correlations through the squared coherence statistic and Lund test discussed in Section 4.1. For example, Figure 11 shows the plots for HSWST (housing start west) variable. Here, periodic correlations should manifest themselves at the index T/s = 80/4 = 20, which is seen more clearly from the Lund test statistic plot. The scree plot is presented in the left panel of Figure 12, but it is not very informative to decide on a number of factors. The information criterion with the penalty function g 1 in (2.20) suggested two factors, which we use here. The other penalty functions g 2 and g 4 suggested only one factor and g 3 suggested four factors. Based on two factors, the estimated factors and loadings are depicted in Figures 12 and 13. Similar overall pattern of loadings amongst quarters suggests DFM or PDFM with common loadings.
The best model is DFM-sPVAR(1;6) showing the smallest h-step-ahead forecasting errors for all h = 1, 2 and 4 (see Table 7). It is remarkable that it only requires 6 non-zero parameters -just 1.37% of a univariate PAR(1) model while giving the smaller h-step-ahead MSE. The estimated sPVAR(1;6) parameters are presented in Figure 14. It shows that most of periodic nature is due to lag-1 dependence of the first factor, and that the second factor makes a significant effect cross-sectionally for the third quarter. The sum of 12-step-ahead forecasting errors of block-DFM as in (4.2) is 838.31 while that of DFMS-sPVAR is 697.56, so PDFM performs better than the block-DFM approach.

Conclusions
In this work, we introduced a periodic dynamic factor model (PDFM) that generalizes the usual dynamic factor model (DFM) by allowing its loadings to depend on "season," and its factors to follow a periodic VAR model with    coefficients depending on "season" (the so-called PVAR model). Two estimation approaches for PDFM were considered: the EM estimation for the so-called exact formulation and the PCA estimation for the approximate formulation. Some asymptotic results were proved for the approximate formulation, by adapting the work on DFM by Doz et al. (2011Doz et al. ( , 2012. The simulation study showed satisfactory performance of the proposed estimation procedures and two real data sets were considered to illustrate the use of PDFMs and related models. A number of open questions related to PDFMs and the standard DFMs were already raised in the text, including extensions to long-memory and nonlinear models, and to weak factors (Remark 2.2), and identifiability issues related to the varying number of factors (Section 2.5). Another interesting direction would be to develop detection and testing methods for the various special cases of PDFMs, as in (1.6). In this work, these various special cases were assessed just based on the out-of-sample forecasting performance.
Therefore, the Fourier-PVAR model can be written in a matrix form as By incorporating possible correlations of ζ n through an estimated covariance matrix Σ of ζ n , the adaptive lasso estimator of Fourier-PVAR coefficients is defined iteratively by where α j is the jth element of α and the weights are proportional to the inverse of coefficients, namely, The covariance matrix at the th iteration is obtained by We have used the 10-fold cross-validation rule to find the tuning penalty parameter λ and the initial estimator of the correlation matrix is given by the identity matrix of dimension q.