Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
59,276 characters · 12 sections · 70 citation commands
Approximate State Space Modelling of Unobserved Fractional Components
\thispagestyle{empty} \setcounter{page}{0}
\paragraph{\bf Abstract.}
We propose convenient inferential methods for potentially nonstationary multivariate unobserved components models with fractional integration and cointegration. Based on finite-order ARMA approximations in the state space representation, maximum likelihood estimation can make use of the EM algorithm and related techniques. The approximation outperforms the frequently used autoregressive or moving average truncation, both in terms of computational costs and with respect to approximation quality. Monte Carlo simulations reveal good estimation properties of the proposed methods for processes of different complexity and dimension.
\paragraph{\bf Keywords.}
Long memory, fractional cointegration, state space, unobserved components.
\paragraph{\bf JEL-Classification.}
C32, C51, C53, C58.
Fractionally integrated time series models have gained significant interest in recent decades. In possibly nonstationary multivariate setups, which arguably bear most potential e.g.\ for assessing macroeconomic linkages, and which are essential for the joint modelling of financial processes, several parametric models have been explored. Among the most popular are the fractionally integrated VAR model Nie2004b, the triangular fractional cointegration model of RobHua2003 and the cointegrated VAR$_{d,b}$ model of Joh2008.
Meanwhile, also models with unobserved fractional components have proven useful, as empirical and methodological work by RayTsa2000, Mor2004, CheHur2006, Mor2007 and LucVer2015 documents. The unobserved fractional components model allows for a generalization of the classic trend-cycle decomposition, where the long-run component is typically assumed to be I(1). As HarTscWeb2019 show, the model can be used to test the I(1) assumption against a fractional alternative. Furthermore, unobserved fractional components allow the formulation of parsimonious models, like factor models, in an interpretable way.
These methods offer a variety of potential applications to empirical researchers. Long-run components of GDP, (un-)employment, and inflation are typically estimated via unobserved components models that restrict the integration order to unity MorNelZi2003, DomGom2006, KliWeb2016. Since there is comprehensive evidence for long memory in these variables HasWol1995, DijFraPa2002, TscWebWe2013 unobserved fractional components may provide new insights regarding the form and persistence of the long-run components. On the other hand, fractionally integrated factor models are constructed straightforwardly using fractional unobserved components. They can be used to assess fractional cointegration relations and for forecasting, as HarWei2018a demonstrate.
Inferential methods for such unobserved fractional components are the subject of this paper. So far, the bulk of empirical work in this field has been conducted in a semiparametric setting, which may be explained by the high computational and implementation cost of state-of-the-art parametric approaches such as simulated maximum likelihood MesKooOo2016. Especially for models of relatively high dimensions or with a rich dynamic structure, there is a lack of feasible estimation methods. Furthermore, in most empirical applications, methods are required to smoothly handle nonstationary cases alongside stationary ones.
We consider a computationally straightforward parametric treatment of fractional unobserved components models in state space form. An approximation of potentially nonstationary fractionally integrated series using finite-order ARMA structures is suggested. This procedure outperforms the more commonly used truncation of fractional processes ChaPal1998 by providing a substantial reduction of the state dimension and hence of computational costs for a desired approximation quality. We derive both, the log likelihood and an analytical expression for the corresponding score. Hence, parameter estimation by means of the EM algorithm and gradient-based optimization make the approach feasible even for high dimensional datasets. In Monte Carlo simulations we study the performance of the proposed methods and quantify the accuracy of our state space approximation. For fractionally integrated and cointegrated processes of different dimensions, we find promising finite-sample estimation properties also in comparison to alternative techniques, namely the exact local Whittle estimator, narrow band least squares and exact state space methods. By using a parameter-driven state space approach, our setup inherits several additional favorable properties: Missing values are treated seamlessly, several types of structural time series components such as trends, seasons and noise can be added without effort, and a wide variety of possibly nonlinear or non-Gaussian observation schemes may be straightforwardly implemented; see Har1991,DurKoo2012.
In this paper we apply the proposed approximation scheme to a $p$-dimensional observed time series $y_t$, which is driven by a fractional components (FC) process as defined by HarWei2018a,
Here, $\varLambda$ is a $p \times s$ coefficient matrix with full column rank, the latent process $x_t = (x_{1t}, ..., x_{st})'$ holds the purely fractional components which are driven by a noise process $\xi_{t}=(\xi_{1t}, ..., \xi_{st})'\sim \mathrm{NID}(0, I)$, and $u_t$ holds the short memory components.
More precisely, while the stationary series $u_t$ is only required to have a finite state space representation, the components of the $s$-dimensional $x_t$ are fractionally integrated noise according to
where for a generic scalar $d$, the fractional difference operator is defined by
and $L$ denotes the lag or backshift operator, $Lx_t = x_{t-1}$. We adapt a nonstationary type II solution of these processes Rob2005 and hence treat $d_j\geq0.5$ alongside the asymptotically stationary case $d_j<0.5$ in a continuous setup.
The fractional unobserved components framework captures univariate and multivariate processes with both long-run and short-run dynamics, fractional cointegration and polynomial cointegration, as well as possibly high-dimensional processes with factor structure. It allows for an intuitive additive separation of long run and short run components, i.e. cyclical and trend components in business cycle analysis, while obtaining similar flexibility as (cointegrated) multivariate ARFIMA models. See HarWei2018a for the relation of the FC model to several other fractional integration setups.
The paper is organized as follows. Section (ref) discusses the state space form, while section (ref) outlines maximum likelihood estimation. In section (ref), the estimation properties are investigated by means of Monte Carlo experiments before section (ref) concludes.
Unlike the stationary long-memory processes considered in the literature, e.g., by ChaPal1998, HsuRayBr1998, HsuBre2003, Bro2007, MesKooOo2016 as well as GraMag2012, our nonstationary type II specification of fractional integration is straightforwardly represented in its exact state space form by setting starting values of the latent fractional process to zero, $x_{jt}=0$ for $t\leq 0$. The solution for $x_{jt}$ is based on the truncated operator $\Delta_+^{-d_j}$ Joh2008 and given by \[ x_{jt} = \Delta_+^{-d_j} \xi_{jt} = \sum_{i=0}^{t-1} \pi_i(-d_j) \xi_{j,t-i}, \qquad j=1,\ldots,s. \] For a given sample size $n$, $x_t$ has an autoregressive structure with coefficient matrices $\varPi_j^{d}$ $= \operatorname{diag}(\pi_j(d_1), \ldots, \pi_j(d_s))$, $j=1,\ldots,n$. Thus, a Markovian state vector embodying $x_{t}$ has to include $n-1$ lags of $x_{t}$ and is initialized deterministically with $x_{-n+1}=\ldots=x_{0}=0$. In principle, this exact state space form can be used to compute the Kalman filter, to evaluate the likelihood and to estimate the unknown model parameters by nonlinear optimization routines. Since the state vector is at least of dimension $s\cdot n$, this can become computationally very costly, particularly in large samples and for a large number $s$ of fractional components, which makes a treatment of the system in its exact state space representation practically infeasible for a wide range of relevant applications.
For non-negative integration orders, note that the $x_t$ can be generalized to have non-zero starting values $x_0$. In that case, $x_t$ is initialized via a diffuse initial state vector. Details on the initialization of nonstationary components are given in Koo1997. Deterministic components in $x_t$ are handled straightforwardly by defining $x_{jt}^* := x_{jt} + \mu_{jt}$ where $\operatorname{Var}(\mu_{jt})=0$ $\forall j = 1,..., s$. Depending on the integration order $d_j$ the contribution of $\mu_{jt}$ on $y_t$ converges to a trend of degree $\lfloor d_j \rfloor$ as $t \rightarrow \infty$.
The literature on stationary long-memory processes has considered approximations based on a truncation of the autoregressive representation, considering only $m$ lags of $x_t$ for $m<n$ in the transition equation (i.e., setting all autoregressive coefficients to zero for $j>m$). Alternatively, the moving average representation has been truncated to arrive at a feasible state space model; see Pal2007, sections 4.2 and 4.3.
Instead, we will apply ARMA approximations to the fractional state vectors, which provide a better approximation quality than the autoregressive or moving average truncation. An ARMA approximation of long-memory processes has been considered in the importance sampling frameworks of HsuBre2003 and MesKooOo2016, but, arguably due to their computational burdens, did not find usage in applied research so far. In our setup, where fractional integration appears in the form of purely fractional components rather than ARFIMA processes, this approach is particularly convenient. In contrast to recent attempts to approximate ARFIMA processes by ARMA ones BasChaPa2001, we do not freely estimate all ARMA parameters but only $d$, and thus retain the original parsimonious parameterization of the process.
As a (nonstationary) approximation of a generic univariate $x_t = \Delta_+^{-d}\xi_t$, we consider the process
for finite $v$ and $w$, where $\varphi:=(a_1,\ldots,a_v,m_1,\ldots,m_w)'$ and all $a_i$ and $m_j$ are made functionally dependent on $d$ to approximate $x_t$ by $\tilde{x}_t$. In order to determine the parameters $\varphi$, we minimize the distance between $x_t$ and $\tilde x_t$, using the mean squared error (MSE) over $t=1,\ldots,n$ as the distance measure. For given $t$, $d$ and $\varphi$, we observe
Hence, the MSE for period $t$ is given by \[ E[ (\tilde{x}_t -x_t)^2] = Var(\xi_t) \sum_{j=0}^{t-1} (\tilde{\psi}_j(\varphi)-\psi_j(d))^2, \] while averaging over all periods for a given sample size $n$ and ignoring the constant variance term yields the objective function for a given $d$,
The approximating ARMA coefficients are thus given by
To obtain the approximating ARMA coefficients in practice, we conduct the optimization (ref) over a reasonable range of $d$, such as $d \in [-0.5; 2]$, for a given $n$. Computational details of the optimization are given in Appendix (ref). Interestingly, for $d < 1$, stationary ARMA coefficients provide the minimum MSE, while for $d\geq1$ we impose an appropriate number of unit roots to enhance the approximation quality.
To illustrate the results we plot the approximating ARMA(2,2) parameters as a function of $d$ for $n=500$; see figure (ref). A closer look at the coefficients reveals that for $d>0$ typically both the autoregressive and the moving average polynomial have roots close to unity which nearly cancel out. For example, to approximate a process with $d=0.75$ we have $(1-1.932L +0.932L^2)\tilde{x}_t = (1-1.285L+ 0.306L^2) \xi_t$, which can be factorized as $(1- 0.999L) (1-0.933L)\tilde{x}_t = (1-0.970L)(1- 0.316L) \xi_t$. Despite their similarity, AR and MA roots do not cancel out for non-integer $d$, since the approximation quality is improved by additional free parameters. For integer integration orders the optimization yields $(1-0.953L)(1-0.477L)\tilde{x}_t = (1-0.953L)(1-0.477L)\xi_t$ for $d=0$ and $(1-L)(1-0.980L)\tilde{x}_t = (1-0.980L)(1+0.001L)\xi_t$ for $d=1$. Consequently, our ARMA-approximation is consistent with the finite representation of inter-integrated processes.
To compare the ARMA($v$,$w$) approximations with $v=w \in \{1,2,3,4\}$ to a truncated AR($m$) process, we contrast the approximating impulse response function $\tilde{\psi}_j$ to the true one, $\psi_j(d)$, for a given $d$. The autoregressive truncation lag $m=50$ is used for our comparison, since this is among the largest values which we consider as feasible in a typical multivariate application. The result of this comparison is shown in figure (ref) for $n=500$ and $d=0.75$. The autoregressive truncation approach gives the exact impulse responses for horizons $j\leq50$, but then tapers off too fast. The ARMA approximations improves significantly over the autoregressive truncation whenever $v=w\geq2$. For orders 3 or 4, the approximation error is even hardly visible. For the moving average truncation, the impulse responses equal zero for horizons exceeding the truncation lag (not shown).
To perform the comparison for different $d$, we plot the square root of the MSE (ref) as a function of $d$ for different approximation methods. For negative integration orders, as shown in figure (ref), the moving average approach clearly outperforms the autoregression, while the ARMA method with orders $v=w>2$ are better. The moving average approximation becomes inaccurate, however, for the case $d>0$, and worse even than the autoregressive method as can be seen in figure (ref). In contrast, the ARMA(3,3) and ARMA(4,4) approximations are well-suited to mimic fractional processes over the whole range of $d$. Further evidence in favor of the ARMA approximation will be presented in the Monte Carlo simulation of section (ref).
Based on these methods we introduce the state space form of the multivariate model (ref), where each $x_{jt}$ is approximated by the ARMA approach. In the following we drop the tilde for the approximation of $x_{jt}$ for notational convenience. To cover the very general case, we allow for residual auto- and cross-correlation by modelling the latent $p$-dimensional short memory process $u_t$ via a stationary state space model, which can capture vector autoregressive, vector ARMA or factor models, among others, and include an additional noise term $\varepsilon_t $. The model can be written in state space form as
where the states may be partitioned into $\alpha'_t =(\alpha^{(1)'}_t,\alpha^{(2)'}_t)$, the states related to the fractional and the stationary components, respectively.
Regarding the fractional part, we define $A_j^d:=\operatorname{diag}(\hat{a}_j(d_1),\ldots,\hat{a}_j(d_s))$ and $M_j^d:=\operatorname{diag}(\hat{m}_j(d_1),\ldots,\hat{m}_j(d_s))$ which contain the approximating AR and MA coefficients of the fractional noise introduced in section (ref), while $A_{j}^d = 0$ for $j>v$ and $M_{j}^d = 0$ for $j>w$. Then $(I - A_1^d L - ... - A_v^dL^v)(I + M_1^d L + ... + M_w^dL^w)^{-1}\tilde{x}_t = \xi_t$. For a minimal state space representation, define $\mu_t := (I + M_1^d L + ... + M_w^dL^w)^{-1}\tilde{x}_t$ such that $(I - A_1^d L - ... - A_v^dL^v) \mu_t = \xi_t$. For $u=\max(v,w+1)$, the first part of the state vector is a $(us)$-dimensional process $\alpha^{(1)'}_t=(\mu_t',\ldots, \mu_{t-u+1}')$. Thus, $\alpha^{(1)}_{t+1} = T^{(1,1)} \alpha^{(1)}_{t} + R^{(1)} \eta_{t}^{(1)}$ with $\eta_{t}^{(1)} = \xi_t$, $R^{(1)'} = (I, 0, \ldots, 0)^{'}$, $Q^{(1,1)} = I$ and \[ T^{(1,1)} =
. \] The observation equation for the fractional part is $\tilde{x}_{t} = \mu_t + M_1^d\mu_{t-1}+ \ldots +M_{u-1}^d \mu_{t-u}$, which enters the observed process $y_t$ through $\varLambda \tilde{x}_{t} = Z^{(1)}\alpha^{(1)}_t$. Thus, the observation matrix for the fractional part is \[ Z^{(1)} =
. \]
For the nonfractional part, we allow a general specification with $\alpha^{(2)}_{t+1} = T^{(2,2)} \alpha^{(2)}_{t} + R^{(2)} \eta^{(2)}_t$ and $u_t = Z^{(2)} \alpha_{t}^{(2)}$, where the distribution of unknown parameters over $T^{(2,2)}$ and $Z^{(2)}$ reflect the choice of the specific model. Without loss of generality, we set $Q^{(2,2)} = \operatorname{Var}(\eta^{(2)}_t) = I$, so that scales and cross correlations of $u_t$ are determined by $Z^{(2)}$. The full state space model (ref) is given by an obvious definition of the system matrices as $Z = (Z^{(1)},Z^{(2)})$, $R' = (R^{(1)'}, R^{(2)'})$, $T = \operatorname{diag}( T^{(1,1)}, T^{(2,2)})$ and $Q = I$. The dynamics are complemented by the initial conditions for the states. From the definition of our type II fractional process we set fixed starting values such as $\alpha^{(1)}_0 =0$, while $\alpha^{(2)}_t$ is initialized by its stationary distribution.
The fractional components $\tilde{x}_t$ do not explicitly appear as states in this representation. However, filtered and smoothed states can be constructed using the relation $\tilde{x}_t = \mu_{t}+\sum_{j=1}^{w} M^d_j\mu_{t-j}$. To obtain conditional covariance matrices for $\tilde{x}_t$, it is more convenient to use an alternative state space form of the ARMA process, where the MA coefficients appear in $R^{(1,1)}$ rather than in $Z^{(1)}$; see DurKoo2012, section 3.4. The current setup, however, is appropriate for estimating the parameters via the EM algorithm which is discussed in the next section.
The EM algorithm was proposed for maximum likelihood estimation of state space models by ShuSto1982 and WatEng1983. Especially in the context of high-dimensional dynamic factor models with possibly more than hundred observable variables, i.e. $p>100$, this method has been found very useful in finding maxima of high-dimensional likelihood functions; see, e.g., QuaSar1993, DozGiaRe2012 and JunKoo2014. After rapidly locating an approximate optimum, the final steps until convergence are typically slow for the EM algorithm, and hence it has been suggested to switch to gradient-based methods with analytical expressions for the likelihood score at a certain step.
We will present these algorithms for our fractional model and thereby extend existing treatments in the literature. For the model represented by (ref), the matrices $T$ and $Z$ both nonlinearly depend on $d$ and other unknown parameters, so that there are nonlinear cross-equation restrictions linking the transition and the observation equation of the system.
The EM algorithm in general consists of two steps, which are repeated until convergence. In the E-step the expected complete data likelihood is computed, where the expectation is evaluated for a given set of parameters $\theta_{\{j\}}$, while the M-step maximizes this function to arrive at the parameters used in the next E-step, $\theta_{\{j+1\}}$. Thus, we define $Q(\theta,\tilde{\theta}) := \operatorname{E}_{\tilde{\theta}}\left[ l(\theta)\right]$, where in this section all expectation operators are understood as conditional on the data $y_1,\ldots,y_n$. In the course of the EM algorithm, after choosing suitable starting values $\theta_{\{1\}}$, the optimization $\theta_{\{j+1\}} = \arg \max_{\theta} Q(\theta, \theta_{\{j\}})$ is iterated for $j=1,2,\ldots$ until convergence.
To state the algorithm for the model defined by (ref) and specified further in section (ref), we follow WuPaiHo1996 to obtain the expected complete data likelihood as
where in our case $Q=I$, while $T$, $Z$ and $H$ are functions of the vector of unknown parameters $\theta$ and a possible dependence of the initial conditions for $\alpha_0$ on $\theta$ has been discarded for simplicity. The conditional moment matrices $A_{\{j\}}$, $B_{\{j\}}$, \ldots, are given in appendix (ref) and can be computed by a single run of a state smoothing algorithm DurKoo2012 based on the system determined by $\theta_{\{j\}}$.
Rather than carrying out the full maximization of $Q(\theta, \theta_{\{j\}})$ at each step, we obtain a computationally simpler modified algorithm. To this end, we partition the vector of unknown parameters as $\theta'=(\theta^{(1)'},\theta^{(2)'})$ where $\theta^{(1)'} = (d',\lambda',\varphi')$, $\lambda$ contains the unknown elements in $\varLambda$, $\varphi$ holds the unobserved parameters for $u_t$ in $T^{(2,2)}$ and $Z^{(2)}$, while the noise variance parameters in $H$ are collected in $\theta^{(2)}$. First, the expectation / conditional maximization (ECM) algorithm described by MenRub1993 in our setup amounts to a conditional optimization over $\theta^{(1)}$ for given variance parameters $\theta^{(2)}_{\{j\}}$ and optimization over $\theta^{(2)}$ for given $\theta^{(1)}_{\{j\}}$. Second, as suggested by WatEng1983, the optimization over $\theta^{(1)}$ is not finalized for each $j$, but rather a single Newton step is implemented for each iteration of the procedure. Neither of these departures from the basic EM algorithm hinders reasonable convergence properties.
A Newton step in the estimation of $\theta^{(1)}$ for given $\theta^{(2)}_{\{j\}}$ yields the estimate in the $(j+1)$-th step
The derivation of (ref) and expressions for $\varXi_{\{j\}}$, $\xi_{\{j\}}$, $g_{\{j\}}$ and $G_{\{j\}}$ can be found in appendix (ref). Finally, the free variance parameters of $H$, collected in $\theta^{(2)}$, are estimated using the derivative of $Q(\theta,\theta_{\{j\}})$ with respect to $H$; see (ref). The estimate is given by the corresponding elements of \[ \frac{1}{n}L_{\{j\}} := \frac{1}{n} \operatorname{E}_{\theta_{\{j\}}}\sum_{t=1}^n \varepsilon_t \varepsilon_t' = \frac{1}{n}(D_{\{j\}} - Z E_{\{j\}}' - E_{\{j\}} Z' + Z F_{\{j\}} Z'). \]
For using gradient-based methods in later steps of the maximization, the likelihood score can be obtained with only one run of a state smoothing algorithm. This has been shown by KooShe1992, who draw on the result \[ \left.\frac{\partial Q(\theta,\theta_{\{j\}})}{\partial \theta}\right|_{\theta_{\{j\}}} = \left.\frac{\partial l(\theta)}{\partial \theta}\right|_{\theta_{\{j\}}}, \] where $l(\theta)$ denotes the Gaussian log-likelihood of the model. Evaluation of the score for our model can therefore be based on (ref) and (ref).
An estimate of the covariance matrix can be computed using an analytical expression for the information matrix. Denoting by $v_t$ and $F_t$ the model residuals and forecast error variances obtained from the Kalman filter, the $i$-th element of the gradient vector for observation $t$ is given by
while the $ij$-th element of the information matrix ${\cal I}(\theta)$ is
see Har1991. To obtain a feasible estimator $\hat{{\cal I}}(\hat{\theta})$, either the expectation term in (ref) is omitted, as suggested by Har1991, or the techniques of CavShu1996 may be used to compute the exact Fisher information. An estimate of the covariance matrix of the estimator is then given by
or by the sandwich form
which is robust to certain violations of the model assumptions; see Whi1982.
The asymptotic theory for maximum likelihood estimation in the fractionally cointegrated state space setup with integration orders $d \in [0, 1.5)$ is derived in HarTscWeb2019 for an exact representation of (ref). As shown there, the approximation error of the Kalman filter that results from ARMA approximations can be calculated via (ref) and is $\operatorname{E}_\theta(\tilde{x}_t - x_t | y_1,...,y_{t-1})$. Hence, it is measurable given the $\sigma$-field generated by $y_1,...,y_{t-1}$, such that an approximation-corrected estimator can be constructed HarTscWeb2019. For this estimator, consistency and asymptotic (mixed) normality is shown. While the approximation-corrected maximum likelihood estimator is computationally feasible for long time series (i.e.\ $n$ large), it is limited to low-dimensional $y_t$. Thus, especially for models where $y_t$ typically holds a large number of observable variables, e.g.\ factor models, the proposed ARMA approximations provide a computationally feasible parametrization of state space models. We compare the performance of the maximum likelihood estimator for our approximate state space model with the approximation-corrected maximum likelihood estimator of HarTscWeb2019 in a Monte Carlo study in section (ref), where it will become clear that the mean squared error of the approximate estimator for $d$ converges to the mean squared error of the exact estimator as $n$ increases.
Our estimation approach can be straightforwardly generalized to additional situations of great practical relevance. To include a treatment of further components causing nonstationarity such as deterministic trends or exogenous regressors, one can use diffuse initialization of one or more of the states which may be based on Koo1997. While we have discussed maximum likelihood estimation under a setting where all data in $y_t$ are available, our algorithms can be generalized for arbitrary patterns of missing data using the approach of BanMod2012. For very high-dimensional datasets, the computational refinements of JunKoo2014 may be used. For common trends of similar persistence, nonparametric averaging methods may turn out to be useful ErgRod2016.
We study the performance of the described methods for a number of stylized processes which are nested in the general setup (ref). The simulation study is designed to answer several questions. Firstly, we assess whether the finite-order ARMA approximation of the state space system performs well as compared to other parametric or semiparametric approaches. Secondly, we assess the feasibility of joint estimation of memory parameters and cointegration vectors in bivariate fractional systems with and without polynomial cointegration, again considering popular semiparametric approaches as benchmarks. Thirdly, the precision of cointegration estimators is studied in case of several cointegration relations of different strengths and for higher dimensions of the observed time series.
For each specification, we simulate $R=1000$ replications and estimate the models using semiparametric estimates for $d$ from the exact local Whittle estimator as starting values for maximum likelihood estimation. The coefficients of the unobserved components can be recovered via the variance of the fractionally differenced observables, since the disturbance terms are standardized. The precision of the estimators is assessed by the root mean squared error (RMSE) criterion or the bias or median errors of the parameter estimators, of state estimates or of out-of-sample forecasts. We vary over different sample sizes $n \in \{250, 500, 1000\}$ which cover relevant situations in macroeconomics and finance.
As the simplest stylized setup of our model, we first assess the fractional integration plus noise case, which has been studied in a stationary setup, e.g., by GraMag2012. For mutually independent $\xi_t$ and $\varepsilon_s$, the data generating process is given by
The fractional integration plus noise model is a special case of (ref) where $\varLambda = \sqrt{q}$, $u_t = \varepsilon_t$, $\mathrm{Var}(\varepsilon_t):=h= 1$, and $\xi_t$, $\varepsilon_t$ are independent. For the signal-to-noise ratio we consider $q \in \{0.5, 1, 2\}$, while the memory parameters $d \in \{ 0.25, 0.5, 0.75\}$ cover cases of asymptotically stationary and nonstationary fractional integration. We estimate the free parameters $d$, $q$ and the noise variance $h$ by maximum likelihood using the state space approach.
We apply different approximations to avoid an otherwise $n$-dimensional state process. Firstly, the ARMA($v$,$w$) approximation given by (ref) and (ref) is considered, setting $v=w \in \{2,3,4\}$. The corresponding estimators are denoted as $\hat{d}_{v,w}$ in the result tables. Secondly, we assess truncations of the autoregressive representation of the fractional process at $m=20$ and $m=50$ lags as suggested in Pal2007, and label these estimators $\hat{d}_{A\!R20}$ and $\hat{d}_{A\!R50}$, respectively. Thirdly, moving average representations as proposed in ChaPal1998 are used, also with a truncation at $m=20$ and $m=50$ lags ($\hat{d}_{M\!A20}$ and $\hat{d}_{M\!A50}$). Furthermore, we employ the exact local Whittle ($\hat{d}_{E\!W}$) estimator of ShiPhi2005 as well as the univariate exact local Whittle approach ($\hat{d}_{U\!E\!W}$) as defined by SunPhi2004, which accounts for additive $I(0)$ perturbations. For both semiparametric estimators of the fractional integration order, we use $m=\lfloor n^{0.65} \rfloor$ Fourier frequencies as a common pragmatic choice. Using other typical values such as $n^j$, $j\in \{0.45, 0.5, 0.55\}$ would not change the results qualitatively, but $n^{0.65}$ is the best choice in most settings considered here. Finally, to grasp the performance of the exact maximum likelihood estimator and to compare our approximate approach with it, we also include the approximation-corrected maximum likelihood estimator of HarTscWeb2019, which corrects for the approximation error induced by ARMA(3,3) approximations.
The root mean squared errors of estimates of $d$ for this setup are shown in table (ref). Not surprisingly, for this stylized process with only three free parameters, the parametric approaches clearly outperform the semiparametric Whittle estimators. For the EW approach, the performance gets worse for more volatile noise processes (lower $q$), which is not the case for the UEW estimator. The bias of the EW estimator is negative due to the additive noise; see table (ref) and also SunPhi2004. In contrast, the UEW estimator is positively biased, independently of $q$. Overall, it has inferior estimation properties, so that we do not show the UEW results for the other data generating processes.
Focusing on the state space approximations, we find that the ARMA approach for $v,w\geq 3$ is always among the best approaches. Overall, the ARMA(3,3) and ARMA(4,4) approximations exert a very similar performance, and their relative performance does not seem to depend on the specification of $d$ and $q$. The truncation methods, in contrast, show mixed results. The moving average approximation tends to dominate the autoregressive one for smaller $d<0.5$, which mirrors the conclusion from GraMag2012 in their stationary setting. However, we find that the autoregression is better whenever nonstationary $d\geq 0.5$ or higher signal-to-noise ratios are considered.
As expected, the exact maximum likelihood estimator of HarTscWeb2019 outperforms the approximation methods for most parameter settings. Considering the computational costs which are about 10 times higher than for the ARMA-approximations with $n = 250$, and about 250 times higher with $n=1000$, the improvements are moderate, however. The median improvement in RMSE over the ARMA(3,3) across parameter setups is 7.7%. As the most extreme scenario, the RMSE can be reduced from 0.132 in the ARMA(3,3) method to 0.102 by the exact estimator for $n =250$, $q=0.5$, $d = 0.25$. As the signal-to-noise ratio $q$ increases, benefits from the approximation-corrected estimator get smaller, and also an increase in the integration orders lowers the benefits from the approximation-corrected estimator. But most interestingly, the RMSE of the ARMA approximations converges to the RMSE of the approximation-corrected estimator as $n$ increases. This indicates that for long time series, where the approximation-correction is particularly costly, it may not even be required.
Directing attention to table (ref) again, we find that the bias for the ARMA approach for $v,w\geq 3$ does not contribute significantly to the estimation errors. Often, it does not appear until the third decimal place. The bias is generally small also for the truncation approaches, but there exist some situations where it is noticeable, mostly for larger $d$. There, larger sample sizes even tend to increase the bias, while higher truncation lags do not always lessen the problem.
We investigate if the results carry over to estimation precision of the fractional components and to forecasting performance of the different approaches. This seems to be the case as table (ref) shows. We restrict attention to the medium signal-to-noise case $q=1$ and apply only the state space approaches. In the upper panel, the fractional component $x_t$ is estimated by a Kalman smoother and the RMSE averages across all in-sample observations and iterations. We find rather small differences between the approaches, especially for small $d$, while ARMA(4,4) and ARMA(3,3) dominate the approximation-based methods in each constellation. The exact method is only slightly superior. The same holds for the forecasting performance, where 1- to 20-step ahead forecasts are evaluated against realized trajectories, again by their RMSE averaging both across horizons and iterations. The differences between the approaches are only slightly more pronounced than above, and again, ARMA(3,3) and ARMA(4,4) are very close to the best-performing exact method. Interestingly, with the low signal-to-noise ratio (not shown in the table), each approach does a poorer job to recover the underlying fractional component, and neither is able to appropriately separate the fractional from the noise component in any case.
In sum, we find good performance of the ARMA approximations. The ARMA(3,3) approach appears sufficient in typical empirical applications. This finding is very appreciable in light of the great reduction in computational effort: A fractional component is represented by 4 states, rather than by 50 in a truncation setup with inferior performance, while an approximation-corrected approach has higher computational costs especially for $n=1000$ even in this very simple setup. Both these alternatives can easily become impractical in more complex situations.
Overall, the differences between the approximations account for a small fraction of the overall estimation uncertainty, even in this stylized setting with high overall estimation precision. Also the benefits of the approximation-corrected approach are limited. Together with the finding of accurate ARMA approximations in section (ref), this suggests that the need of approximations might not be a serious obstacle to the state space modelling of fractional unobserved components.
The performance of the state space approach in estimating fractionally cointegrated systems is studied in a bivariate process with short-run dynamics,
with $\Lambda_{11}=\Lambda_{21}=1$, $\Gamma_{11}= \Gamma_{22} =c$, $\Gamma_{21} = c \cdot e$, and $\phi_1 = \phi_2 = 0.5$. This implies that the true cointegration vector $B := \varLambda_\perp = (1, -1)'$, $B'
' \sim I(0)$, where the first entry was normalized to one. Again the innovations are mutually independent. Note that $u_{1t}=c z_{1t}$, $u_{2t}=(c\cdot e)z_{1z}+cz_{2t}$, which allows for an interpretation of \eqref{eq:DGP2} as a fractionally cointegrated setup with cross- and autocorrelated short-run dynamics. We vary over values of the fractional integration order $d \in \{ 0.25, 0.5, 0.75\}$. The perturbation parameter $c\in\{0.5,1,2\}$ controls the signal-to-noise ratio and short-memory correlation between the processes is introduced, which will be governed by different values of $e \in \{0, 0.5, 1\}$. Cases where $\Gamma_{11} \neq \Gamma_{21}$ could be considered straightforwardly.
Here and henceforth, we apply the ARMA(3,3) approximation for maximum likelihood estimation of the unknown model parameters. In the current setup, the latter consist of the eight entries in $\theta' = (d, \phi_1, \phi_2, \varLambda_{11}, \varLambda_{21}, \varGamma_{11},\varGamma_{21},\varGamma_{22})$, where $\varGamma_{ij}$ is the loading of $z_{jt}$ on $y_{it}$, while the variance parameters are normalized to achieve identification. Starting values for the AR parameters are obtained by fitting an autoregressive model for the difference $y_{1t}-y_{2t}$. To contrast the properties to standard semiparametric approaches again, we apply the EW estimator componentwise to the univariate processes and investigate the mean of the univariate estimates. For the cointegration relation we apply the narrow-band least squares estimator which has been studied by RobMar2001 in the nonstationary single equation case and by Hua2009 in a setup with cointegration subspaces HuaRob2010, HarWei2018a. We follow the literature which suggests to use a small number of frequencies and choose $\lfloor n^{0.3} \rfloor$, amounting to 5, 6 and 7 frequencies for our sample sizes.
Since the cointegration vectors are not identified without further restrictions, we investigate the angle $\vartheta$ between true and estimated cointegration spaces. Nie2010 provides an expression for the sine of this angle, which is given in our framework by
where $\hat{B}$ is an estimated cointegration matrix and $\| A \|$ is the Euclidean norm of $A$. In the current bivariate setup with one cointegration relation, we have $\hat{B} = \hat{\varLambda}_{\bot}$ for the maximum likelihood estimator and $\hat{B}_{N\!B} = (1, -\hat{\beta}_{N\!B})'$ for the narrow-band least squares estimator $\hat{\beta}_{N\!B}$ applied to $y_{1t} = \beta y_{2t} + \text{error}$. Values of $\sin(\vartheta)$ closer to zero indicate preciser estimates and thus we compute the corresponding root mean squared error criterion as the square root of $\frac{1}{R}\sum_{i=1}^R \sin(\vartheta^i)^2$ in what follows. To get some intuition for the bivariate case, estimating a true value $B= (1,-1)'$ by $\hat{B} = (1,-1.1)'$ would result in a loss of $\sin(\vartheta) \approx 0.05$.
In table (ref) we show root mean squared errors for memory parameters ($\hat{d}^{M\!L}$ and $\hat{d}^{E\!W}$) and evaluate estimated cointegration spaces (by $\vartheta^{M\!L}$ and $\vartheta^{N\!B}$) applying either the maximum likelihood or the semiparametric technique, respectively. Consider the case $e = 0$ first. Regarding the memory estimators, we find relatively large errors for this data generating process, with root mean squared errors frequently around 0.2 or larger, most prominently when the variances of the short-memory processes are large ($c=2$). The Whittle estimator often performs better than maximum likelihood, especially for smaller $c$ and $d$ and in smaller samples.
For estimating the cointegration space, however, the state space approach appears worthwhile and outperforms narrow band least squares in most constellations. Not surprisingly, strong cointegration relations ($d=0.75$) are precisely estimated, as is cointegration with small short-memory disturbances ($c=0.5$). While the relative merits of maximum likelihood are unchanged for different cointegration strengths, we find that strong perturbations are better captured by the state space estimators. For $c=2$, the RMSE of the semiparametric approach exceeds the parametric RMSE by about $70\%$ in some cases.
Short memory correlation as introduced through $e>0$ overall decreases the precision of the memory estimators. Interestingly, however, the performance of the cointegration estimators improves when $e>0$ is considered. This is the case for both the maximum likelihood and the narrow band approach. To gain some insights into this finding, we assess the typical signed errors of the cointegration estimates. To this end, we consider a normalization of the cointegration vectors as $(1, -\beta)$, and assess estimated $\beta$ for both approaches. Note that the narrow-band least squares estimator estimates $\beta_{NB}$ directly, whereas $\beta_{ML}$ is computed via $\beta_ {ML}=-\hat{\Lambda}_{21}/\hat{\Lambda}_{11}$. For $\hat{\Lambda}_{11}$ small, the estimator becomes imprecise. Therefore, it is informative to compute an outlier-robust measure of the typical signed deviation. The median errors ($\operatorname{median}_i(\hat{\beta}_j^i)-\beta_j$) for this data generating process are shown in table (ref).
The typical deviations for the narrow band estimates exert a negative median bias of the estimates. A positive correlation between the short-memory components appears to work in the opposite direction so that the negative bias is reduced. In contrast, we find that the maximum likelihood estimators are essentially median-unbiased. Here, correlation between the short-memory components may improve the distinction between short and long-memory components and hence reduce variability.
A further simulation setup extends the model setup (ref) by introducing contemporaneously correlated $\xi_t \sim \mathrm{NID}(0, S)$, and also allowing for polynomial cointegration through perfectly correlated $\xi_t$. Polynomial cointegration refers to a situation where lagged observations nontrivially enter a cointegration relation; see GraLee1989 as well as Joh2008 for nonfractional and fractional treatments, respectively. To motivate polynomial cointegration in terms of our model, assume for simplicity $d_1 > d_2 > ... > d_s$. Let $\Lambda^{(1:(m-1))}$ hold the first $m-1$ columns of $\varLambda$, and let $\varLambda_\perp^{(1:(m-1))}$ be its orthogonal complement. Then $\varLambda_\perp^{(1:(m-1))'}y_t \sim I(d_{m})$ annihilates the first $m-1$ common unobserved components $x_{1t}, ..., x_{m-1,t}$. If a vector $\gamma$ exists, such that $\gamma' (y_t' \varLambda_\perp^{(1:(m-1))} ,\ \Delta^b y_{it})'$ is integrated of a lower order than $d_{m}$ for any $b$ and any $y_{it} \sim I(d_k)$, $k \in \{1,...,m-1\}$, then polynomial cointegration occurs. Whenever $| \mathrm{Cor}(\xi_{kt}, \xi_{mt})| = 1$, also $x_{mt}$ and $\Delta^{d_k-d_m}x_{kt}$ are perfectly correlated, and hence there exists a linear combination $\gamma' ( y_t' \varLambda_\perp^{(1:m)},\ \Delta^{d_k - d_m} y_{it} )'$ with a smaller integration order than $d_{m}$.
We consider
where $\Lambda_{11}=\Lambda_{21}=1$, $\Lambda_{12} = - \Lambda_{22} = a$, $d_1 > d_2$, and where we drop the assumption of orthogonal long-run shocks and allow for $\mathrm{Var}(\xi_t)=Q \neq I$. Correlation between the innovations to the fractional processes is introduced through the parameter $r$. Besides the standard setting $r=0$, we refrain from the assumption of independent components for $r=0.5$, while $r=1$ amounts to $\xi_{1t}=\xi_{2t}$ which is the case of polynomial cointegration since there is a second nontrivial cointegration relation in $(y_{1t}, y_{2t}, \Delta^{d_1-d_2} y_{2t})'$. Combinations of $d_2\in \{ 0.2, 0.4\}$ and $d_1\in \{ 0.6, 0.8\}$ contrast relatively weak and strong cases of cointegration, while the importance of the component $x_{2t}$ varies with $a\in \{0.5,1,2\}$. We treat $\theta = (d_1, d_2, \varLambda_{11},\varLambda_{21},\varLambda_{12}, \varLambda_{22},r, h_{11}, h_{22})'$ as free parameters, but also investigate estimates imposing the singularity $r=1$ when it is appropriate. Starting values for the fractional integration orders are obtained via the exact local Whittle estimator as in the preceding sections, where we consider the sum and the difference of $y_{1t}$ and $y_{2t}$ to estimate $d_1$ and $d_2$. Initial values for $r$ are obtained from the covariance of the fractionally differenced processes $\mathrm{Cov}(\Delta^{d_2}(y_{1t}-y_{2t}), \Delta^{d_1}(y_{1t}+y_{2t}))$.
Consider the results for $r=0.5$ first. The root mean squared errors, shown in table (ref), include estimators of cointegration spaces as above (evaluated by $\vartheta_1^{M\!L}$ and $\vartheta_1^{N\!B}$ in the table). Now, there are two memory parameters to be estimated either by maximum likelihood ($\hat{d}_1^{M\!L}$ and $\hat{d}_2^{M\!L}$) or by the Whittle approach ($\hat{d}_1^{E\!W}$ and $\hat{d}_2^{E\!W}$). Semiparametric estimates of $d_2$ are obtained from the narrow band least squares residuals. The table also contains the maximum likelihood estimate of the correlation parameter $r$ ($\hat{r}^{M\!L}$).
For most parameter settings, we observe that the parametric memory estimators perform satisfactorily. They outperform the semiparametric approach whenever there is strong influence of the $x_{2t}$ components ($a=2$), most pronouncedly in larger samples. Also regarding cointegration estimators, higher values of $a$ favor the parametric method. The correlation parameter is estimated with increasing precision in larger samples, while also the strength of the cointegration relation is relevant for this estimator. For $d_1=d_2$, the correlation parameter (and also certain elements of $\varLambda$) would not be identifiable, and hence setups with small difference $d_1-d_2$ are problematic.
For $r=1$, we additionally consider the properties of estimators for the polynomial cointegration relation. To evaluate estimators of the polynomial cointegration spaces, note that the cointegration space leading to the highest memory reduction in $(y_{1t}, y_{2t}, \Delta^{d_1-d_2} y_{2t})'$ is the orthogonal complement of the span of
where $\varLambda^{(j)}$ refers to the $j$-th column of $\varLambda$. This cointegration subspace is estimated replacing all entries in (ref) by their maximum likelihood estimates, where $r=1$ is imposed. For the narrow band least squares estimator, this space is determined by the span of $(1, -\hat{\beta}_1,-\hat{\beta}_2)'$, where the coefficients are narrow band least squares estimates from $y_{1t} = \beta_1 y_{2t} + \beta_2 \Delta^{d_1-d_2} y_{2t} + \text{error}$ with $d_1$ and $d_2$ replaced by local Whittle estimates. Estimators for this second (polynomial) cointegration relation are evaluated analoguously to (ref) where now (ref) takes the role of $\varLambda$ and the resulting angle is denoted by $\theta_2$.
In table (ref), the corresponding root mean squared errors are given. The elementary cointegration space is estimated by the unrestricted estimator (see $\vartheta_1^{M\!L}$) and the restricted estimator (see $\vartheta_1^{R\!M\!L}$, imposing $r=1$) with a very similar precision. This is in accordance with the notably precise estimation of $r$ in this case. The parametric estimators of both cointegration spaces are again better than semiparametric approaches (1) in large samples and (2) when a strong second fractional component is present. Overall, the results suggest that polynomial fractional cointegration analysis is feasible in our setup, while the maximum likelihood approach has reasonable estimation properties at least for larger sample sizes.
Until now, we have considered one- or two-dimensional processes in our simulations which limits the empirical relevance of the findings so far. We claim that modelling high-dimensional time series constitutes a strength of our approach, at least if suitably sparse parametrizations with factor structures are empirically reasonable. As a second generalisation compared to the previous setups, we consider situations where two or more cointegration relations exist and where these may be of different strength, i.e., where the reduction in memory through cointegration differs among relations. The latter situation has been studied under the label of cointegration subspaces, among others by HuaRob2010 and HarWei2018a.
To assess the performance in this situation, consider the process
where $\Lambda_{i1}=a$, $\Lambda_{i2} = a \cdot (-1)^{i+1}$ $\forall i = 1,...,p$, and with mutually independent noise sequences. We now vary over the dimension $p \in \{3,10,50\}$, while again combinations of $d_1\in \{ 0.2, 0.4\}$ and $d_2\in \{ 0.6, 0.8\}$ are considered. The parameter $a\in \{0.5,1,2\}$ gives the relative importance of the fractional components and hence plays the role of a signal-to-noise ratio. We estimate $d_j$, $\varLambda_{ij}$, $h_i$ for $j=1,2$ and $i=1,\ldots,p$ as free parameters. Starting values for $d_1$ and $d_2$ are obtained as in section (ref).
Along with the memory estimates, we show results for estimating the $p-1$ cointegration relations reducing the memory from $d_1$ to $d_2$ (the first cointegration subspace) which is evaluated by the angle $\vartheta_1$ between $\varLambda^{(1)}$ and the cointegration matrix estimate $\widehat{B}_1$. Additionally, the $p-2$ cointegration relations reducing the memory from $d_1$ to $0$ (the second cointegration subspace) are evaluated by the angle $\vartheta_2$ between $\varLambda$ and $\widehat{B}_2$. The cointegration matrices are straightforwardly obtained for the maximum likelihood approach by the orthogonal complements of $\widehat{\varLambda}^{(1)}$ and $\widehat{\varLambda}$, respectively. The narrow-band least squares method estimates cointegration matrices under specific normalizations as above. Estimating the first subspace, we construct $\hat{B}_1$ to have free entries $-\hat{\beta}_{2}$, \ldots, $-\hat{\beta}_{p}$ in the first row and a $p-1$ identity matrix below, such that $\beta_{j}$ is obtained from $y_{jt} = \beta_j y_{1t} + \text{error}$ for $j=2,\ldots,p$. In the estimation of the second subspace, we have two free rows in $\hat{B}_2$ which are given by $(-\hat{\beta}_{13}$, \ldots, $-\hat{\beta}_{1p})$, and $(-\hat{\beta}_{23}$, \ldots, $-\hat{\beta}_{2p})$, respectively, and can be estimated from $y_{jt} = \beta_{1j} y_{1t} + \beta_{2j} y_{2t} +\text{error}$ for $j=3,\ldots,p$.
In table (ref), results are shown for $a=0.5$ while the other specifications yield qualitatively similar outcomes. The process allows for a precise estimation of both $d_1$ and $d_2$ by maximum likelihood. An increasing dimension $p$ leads to a better estimation by maximum likelihood which is not the case for the Whittle technique. The semiparametric Whittle estimates are obtained by averaging univariate estimates for $d_1$ and using narrow band least squares residuals to estimate $d_2$. Notably, the estimates of $d_2$ hardly improve with larger $n$, which can be explained by a specific shortcoming of the single equation approach: The univariate regression errors may each have integration orders of $d_2$ or lower. In our case, lower orders prevail for $y_{jt} = \beta_j y_{1t} + \text{error}$ with $j$ odd, due to the special structure of $\varLambda$. Knowledge about this specific structure is not exploited by both methods, however, to keep the simulation scenario realistic.
Also regarding the estimation of the cointegration spaces, maximum likelihood is superior. Both parametric and semiparametric estimators have smaller errors for higher dimension, whereas this “blessing of dimensionality” is more pronounced for the state space approach. Generally, the ratio between the maximum likelihood RMSE and the semiparametric RMSE decreases for larger $p$.
Not surprisingly, the case with strongest basic cointegration (large difference $d_1 - d_2$, which implies a great reduction of persistence when $x_1$ is projected out) is the one with highest precision in estimating the first cointegration subspace. For estimating the second subspace, a slightly different logic applies, with a larger $d_2$ supporting the estimation. E.g., in the case $d_1 =0.6$ and $d_2=0.4$ a higher precision is achieved than for $d_1 =0.6$ and $d_2=0.2$. Overall, we find that our approach profits from imposing the factor structure which is not the case for the benchmark methods applied in this comparison.
We have proposed estimation methods for nonstationary unobserved components models which are computationally efficient and provide a good approximation performance. These may be relevant for a wide variety of applications in macroeconomics and finance, as HarWei2018a have illustrated. Further work is needed to assess the performance of the methods in different, possibly very high-dimensional, settings.
The research of this paper has partly been conducted while Roland Weigand was at the University of Regensburg and at the Institute for Employment Research (IAB) in Nuremberg. Very valuable comments by Rolf Tschernig, Enzo Weber, and two anonymous referees, are gratefully acknowledged. Tobias Hartl gratefully acknowledges support through the projects TS283/1-1 and WE4847/4-1 financed by the German Research Foundation (DFG).