Limit Theory under Network Dependence and Nonstationarity

Limit Theory under Network Dependence and Nonstationarity

Introduction

This set of lecture notes is about systems of predictive regressions under network dependence. The first question that we will address is: why should we model network dependence? Our answer to this issue is threefold. First, under both network dependence and time series nonstationarity (using a local-to-unity parametrization) conventional asymptotic theory results are not directly applicable unless we understand the required regularity conditions and modifications that are needed.

Time Series as a Stochastic process

Consider $\left\{ X_{\theta} \right\}_{ \theta \in \Theta }$, a stochastic process with index set $\Theta$. This is a family of random elements defined on a probability space $( \Omega, \mathcal{F}, \mathbb{P} )$. More formally, we have a mapping $X: \Theta \times \Omega \mapsto S$, and thus we can write $X( \theta, \omega )$ so that the element is a function of $\theta$ (possibly time or space) and $\omega$ is randomness. As we shall see in the next section, we can write down a coordinate representation of the probability space so that $\Omega = R_{ \infty } \left( \times_{0}^{\infty} R \right)$ such that $\omega = \left( x_0, x_1,... \right) \in R_{ \infty }$. Here $\omega$ is a typical trajectory or history of the process. We have the following natural maps

The above map picks off the $t'$th element and records a sample path as we change $t$.

The above map produces a random variable that lives on $R$ (the state space) which in other words is a real valued random variable. This perspective is useful because it makes explicit what is the most distinctive feature of time series analysis. The fact that in making inferences about the probability law $\mathbb{P}$ that governs the process the evidence we have available takes the form of a single history $\omega$. In practise, of course, the situation is even worse: we only have a finite stretch of $\omega$ such as $\left( x_0, ..., x_n \right)$ on which to base inference.

Proof. Suppose S is ergodic and take the random variable $1_{G}$. Then, by the ergodic theorem,

The above aspects are considered to be the fundamental assumptions in time series models. However, often these assumptions could vary depending on both the econometric model under consideration as well as on the empirical data application. For example, when modelling stock returns the presence of serial correlation has an important interpretation in the financial economics literature. Specifically, incorporating serial correlation in models of stock prices and expected returns can provide additional evidence on the validity of the efficient market hypothesis and explain other market anomalies.

Further stylized facts of financial returns include the presence of heavy-tailedness and volatility clustering which implies that related econometric assumptions require the existence of fourth moments for the underline distribution function of the model innovations. Further assumptions can include the presence of heteroscedastic errors by imposing a parametric GARCH specification as well as other type of structures such as autoregressive errors. Although the main assumption regarding the innovation sequence of the the stationary AR$(1)$ time series model is that $u_t$ is an i.i.d sequence with mean zero and known variance, this assumption can be further relaxed to incorporate weakly or strongly dependent error sequences. By imposing a linear process representation then we can introduce the notation of weakly dependent errors.

Following the work of granger1978invertibility for a linear stationary process it holds that

for some time series $X_t$ across all $k$.

Therefore, departures from the above expression indicate the presence of nonlinearity in time series. In particular, using the squared residuals from a linear model we can apply a standard Box-Ljung Portmanteau test for serial correlation (see, lee1993testing). In other words, the BLP test is sensitive to departures from linearity in mean.

Martingale Theorems

Covariance Matrix Estimation

Long-Run variance

Denote with $\bar{w} = n^{-1} \sum_{t=1}^n w_t$. When the series is strict white noise (that is, zero mean iid), then the central limit theorem implies that $n^{1/2} \bar{w} \overset{d}{\to} \mathcal{N} (0, V_w)$, where $V_w = \mathsf{Var} (w)$. Furthermore, when the series $w_t$ exhibit serial dependence but is a martingale difference with respect to its own past then this convergence in distribution still holds. However, if the series $w_t$ exhibits serial correlation, the asymptotic normality still holds, but the asymptotic variance becomes

which differs from the variance var$(w)$. Since we have that,

Moreover, we denote the autocovariance of $w_t$ at delay $j$ by $\rho (j) = \mathbb{E} \big[ w_t w_{t-j}^{\top} \big]$. Using that $\rho (-j) = \rho(j)^{\top}$ and the Toeplitz lemma, we obtain

Notice that the above expression for $V_w$ is termed the long-run variance. If instead of zero $\mathbb{E} [ w_t ]$ we consider possibly non-zero $\mathbb{E} \big[ m(w_t) \big]$ for some function $m(.)$, the long-run variance

has the same representation as before, but the autocovariance of $m (w_t)$ at delay $j$ takes the form $\rho (j) := \mathbb{E} \big[ \big( m (w_t) - \mathbb{E} \left[ m (w_t) \right] \big) \big( m (w_{t-j}) - \mathbb{E} \left[ m (w_t) \right] \big)^{\prime} \big]$. Then, the long-run covariance between, say, $m( w_t )$ and $\mu ( w_t )$ is defined similarly such that

is the cross-covariance between $m ( w_t )$ and $\mu ( w_t )$ at delay $j$. Consistent estimation of the long-run variance is not straightforward. Specifically, each autocovariance $\rho(j)$ can be consistently estimated by

for non-negative $j$, and by $\widehat{\rho(j)} = \left( \widehat{\rho (-j)} \right)^{\prime}$ for negative $j$. A natural estimator of $V_m$ is

However, this estimator is not consistent because it is constructed using covariance estimates whose number increases proportionately with the sample size. Therefore, a possible solution is to truncate the summation at $-p$ and $p$, where $p < n - 1$ such that

Due to an ill-posed estimation problem, since the positive-definiteness property of the sample covariance matrix is not guaranteed, NW proposed the following modified version of the covariance matrix:

where the weights ensure that this NW estimator is positive semi-definite. The use of the long-run variance estimators is widespread in estimation and inference for models with dependent data. In particular, $V_n ( \bar{\theta} )$ is a consistent estimator of the long-run variance of the moment function $m ( w_t, \theta_0 )$

Define with

Then, the feasible optimal weighting matrix is then an estimate of $V_m$ which also uses a preliminary estimator $\bar{\theta}$ in place of $\theta_0$. Similarly, for the OLS estimator in a linear time series regression $y_t = x_t^{\prime} \theta_0 + e_t, \mathbb{E} \big[ e_t | x_t \big] = 0$, the asymptotic variance of $\hat{\theta}$ can be estimated using the long-run variance estimator of $x_t \hat{e}_t$, where $\hat{e}_t$ are the OLS residuals of the model. Moreover because, the long-run variance estimators automatically take into account the serial correlation and conditional heteroscedasticity in $e_t$, these estimators are often referred to as HAC estimators (see, newey1986simple, andrews1991heteroskedasticity and andrews1992improved).

Consistent Covariance Matrix Estimation for Linear Processes

Consistency of kernel estimators of the long-run covariance matrix of a linear process is established under weak moment and memory conditions.

We consider the estimation of $\Omega = \displaystyle \underset{ T \to \infty }{ \text{lim} } \frac{1}{T} \sum_{t=1}^T \sum_{s=1}^T \mathbb{E} \left( V_t V_t^{\prime} \right)$, the long-run covariance matrix of $V_t$.

Consistency of kernel estimators of $\Omega$ is established under weak conditions on $\left\{ C_j : j \geq 1 \right\}$ and $\left\{ e_t : t \in \mathbb{Z} \right\}$. We consider a kernel estimator of the form

where $k(.)$ is a measurable kernel function and $\left\{ b_T: T \geq 1 \right\}$ is a sequence of bandwidth parameters. The corresponding estimator of $\Omega$ is given by

which can be written as $\hat{\Gamma}_T + \hat{\Gamma}_T^{\prime} + \hat{\Sigma}_T$, where $\hat{\Sigma}_T = \sum_{ t= 1}^T V_t V_t^{\prime}$. Notice that $\sqrt{n} \left( \hat{ \theta }_n - \theta \right) = \mathcal{O}_p(1)$ implies that $\hat{ \theta }_n \overset{p}{\to} \theta $ with a $\sqrt{n}-$rate of convergence. The order of convergence indicates the rate at which the distance of the sample estimator to the population parameter get smaller as the sample size increases, $n \to \infty$, i.e., $\hat{ \theta }_n = \mathcal{O}_p( \frac{1}{ \sqrt{n} })$. According to kiefer2000simple and kiefer2002heteroskedasticity,

In other words, the asymptotic distribution is a $k-$variate normal distribution with mean 0 and variance-covariance matrix $V = Q^{-1} \Omega Q^{-1}$. Then, the asymptotic distribution of $\hat{\beta}$ can be used to test hypothesis about $\beta$. However, to do this an estimate of $V$ is required. In particular, a natural estimate of $Q$ is

Therefore, $\Omega$ can be estimated by a HAC estimator, $\widehat{\Omega}$. Letting $\hat{u}_t$ be the residuals of the transformed model, the HAC estimate would use $\hat{v}_t = \tilde{F}_t ( \hat{\beta} ) \hat{u}_t$ to estimate nonparametrically the spectral density of $v_t$ at frequency 0, and hence $\Omega$. A typical estimator takes the form

where $k(x)$ is a kernel function satisfying $k(x) = k(-x)$, $k(0) = 1$, $| k(x) | \leq 1$ continuous at $x = 0$ and $\int_{- \infty}^{ \infty } k^2 (x) dx < \infty$. Moreover, the tuning parameter, $s(T)$, is often called the truncation lag or bandwidth. A typical condition for consistency of $\widehat{\Omega}$ is that $s(T) \to \infty$ as $T \to \infty$ but $s(T) / T \to 0$. Therefore, to test hypotheses about $\beta$ in the standard approach, $\widehat{V} = \widehat{Q}^{-1} \widehat{\Omega} \widehat{Q}^{-1}$ is used to transform $\sqrt{T} \left( \hat{\beta} - \beta_0 \right)$ to obtain

In a new approach hong2023kolmogorov (see also sun2022adjusted), use a similar method, but transform $\sqrt{T} \left( \hat{\beta} - \beta_0 \right)$ in such a manner that the asymptotic distribution no longer depends on unknown parameters. The essential difference between the two approaches is that their approach does not require an explicit estimate of $\Omega$ and takes the additional sampling variation associated with not knowing the covariance matrix into account in the asymptotic approximation. On the other hand, HAC estimates treat the variance-covariance matrix as known asymptotically and ignore the impact of finite-sample variability from $\widehat{\Omega}$ on the distribution of test statistics. Thus, a practical limitation of the HAC approach is that to obtain HAC estimates, a truncation lag for a spectral density estimator must be chosen. Although asymptotic theory dictates the rate at which the truncation lag much increase as the sample size grows, no concrete guidance is provided. Consider the following scaled partial sum empirical process:

Their approach provides an elegant solution to this practical problem by avoiding a consistently estimate of the variance-covariance matrix and thus removing the need to choose a truncation lag. In other words, a data-dependent transformation is applied to the NLS estimates of the parameters of interest. This transformation is chosen such that it ensures that the asymptotic distribution of the transformed estimator does not depend on nuisance parameters. The transformed estimator can then be used to construct a test for general hypotheses on the parameters of interest.

Continuous Asymptotics

Before considering the continuous time asymptotics we consider the following example which can be found on page 328 in the book of csorgo1997limit.

Consider a triangular array of random variables $\left\{ \left\{ y_{nt} \right\}_{t=1}^{T_n} \right\}_{n = 1}^{ \infty }$. The triangular array then provides a formal framework within which $h$ may vary and by means of which we may investigate limiting behaviour as $h_n \to 0$. Let $S_{ ni } = \sum_{ j = 1 }^i u_{nj}$ such that $1 \leq i \leq T_n$ with $S_{n0} = 0$.

Furthermore, we form the random function as below

As $n \to \infty, Y_n(r)$ converges weakly to a constant multiple of a standard Wiener process.

We define with

Then, we have the t-statistic given by

where an unbiased estimator for the variance is given by

Furthermore, we have that

and this expansion may be verified directly by developing an expansion for the characteristic function of $\int_0^1 X_T (r)^2 dr$.

The Ornstein–Uhlenbeck process

Hazard rates processes are functionals of piecewise-continuous Gaussian martingales, such as the OU processes and the pure jump process.

\paragraph{Modelling the item state by a stochastic process}

The resulting stochastic differential equation is given by

\paragraph{Modelling item failure rate by a stochastic process}

Dynamic environments exert random stresses on the item. An approach for modelling the lifelength of items in a dynamic environment is to describe the item's failure rate by a stochastic process. Moreover, there is an interesting relation between a hazard rate process and a doubly stochastic Poisson process (known as the Cox process).

Let T denote the lifelength of the item and let $\left\{ \lambda(s), s \leq 0 \right\}$ be a non-negative, real-valued, right-continuous process.Then, $\left\{ \lambda(s), s \leq 0 \right\}$ is said to be the hazard rate process of T if

Consequently,

Higher-Order Autoregressive Processes

AR(p) Time Series Regression

An autoregressive process of order AR$(p)$ is defined by

for any $t \geq 1$ or equivalently, in a compact form we have that $X_t = \theta^{\top} \Phi_{t-1} + \varepsilon_t$, where $\theta = \big( \theta_1,..., \theta_p \big)^{\top}$ is a vector parameter, $\Phi_0$ is an arbitrary initial random vector and $\Phi_t = \big( X_t,..., X_{t-p+1} \big)^{\top}$ and $\varepsilon_t$ is a strong white noise having a finite variance $\sigma^2$. Then, the corresponding characteristic polynomial is defined by

where $E_t = \big( \varepsilon_t, 0, ..., 0 \big)^{\top}$ is a $p-$dimensional noise. Notice that it is well-known that the stability of this $p-$dimensional process is closely related to the eigenvalues of the companion matrix that we will denote and arrange in increasing order as below

The Autoregressive and Moving Average Processes (ARMA)

A stochastic process, or time series, can be repeated as the output resulting from a white noise input, $\varepsilon_t$,

where $X_t$ is considered to be the output series after demeaning e.g., $X_t = \tilde{X}_t - \mu$. In other words, the general linear process representation allows to represent the output of a time series, as a function of the current and previous value of the white noise process, $\varepsilon_t$, which may be represented as a series of shocks. The autocorrelation function of a linear process is given by

Thus, the backward shift operator, $B$, is defined as $B X_t = X_{t-1}$ and $B^j X_t = X_{t-j}$. The autocorrelation generating function may be written as $ \gamma (B) = \sum_{j= - \infty}^{+ \infty } \gamma_k B^k$. Therefore, the strict stationarity assumption holds if and only if the $\phi$ weights of a linear process must satisfy that $\phi(B)$ converges on or lies within the unit circle. Then, an autoregressive, AR, model, the current value of the time series and a random shock $\varepsilon_t$, such that

Therefore, the autoregressive operator of order $P$ is given by

Furthermore, the device proposed byphillips1992asymptotics, is suitable to derive invariance principles and functional central limit theorems in different econometric environments. In particular, gronneberg2019partial focus on establishing limit results for partial sum processes based on ARMAX residuals. The following lemma from gronneberg2019partial adapts the Phillips-Solo device although unlike phillips1992asymptotics, the authors do not impose a time dependence structure, and allow two dimensional arrays.

Therefore, by the Beveridge-Nelson decomposition as stated in Lemma 1 of phillips1992asymptotics,

The second result is stated as below

Denote with $\tilde{Z}_{n,t} := \sum_{j=0}^{\infty} \left( \sum_{k=j+1}^{\infty} c_k \right) Z_{n,t-j}$. Then, we have that

The above result holds due to the following property

Let $\eta > 0$ be given. Then, by definition of supremuem, as $n \to \infty$, we have that

Let $\Gamma(B) = \sum_{j=0}^{\infty} \omega_j B^j$ such that

where $\Gamma(B)\tilde{S}_{n,t}$ is well-defined. Then, addition and subtraction of $\Gamma(1) \sum_{t= 1}^{ \floor{ \lambda n} }\tilde{S}_{n,t}$ with the triangle inequality yields the following result

Moderate Deviations for Multivariate Martingales

According to grama2006asymptotic, limit theorems for probabilities of moderate deviations for sums of independent random variables have been extensively examined before. For martingales, however the number of results on moderate deviations in the literature is still rather limited. Moderate deviation results are related to bounds on the rate of convergence in the central limit theorem. In the multivariate case, quadratic characteristics are matrix valued processes. Therefore, to establish uniform estimates of the rate of convergence in the multivariate martingale CLT, moderate deviations require special handling.

\paragraph{Univariate Case}

Let $\left( \xi_{nk}, \mathcal{F}_{nk} \right)_{ 0 \leq k \leq n}$ be a square integrable martingale difference sequence with $\xi_{n0} = 0$. We set

Denote by $\langle X^n \rangle$ the quadratic characteristic of the martingale $X^n$, such that

\paragraph{Multivariate Case}

Denote by $\langle X^n \rangle$ the quadratic characteristic of the martingale $X^n$, i.e.,

Moreover, it is well known that the distributional properties of the time series $( X_t )_{ t \geq - d + 1 }$ are strongly dependent on the value of $\rho \in \mathbb{R}^d$ are are entirely different for different $\rho-$regions in $\mathbb{R}^d$. One of the basic results states that there exists a stationary distribution of $\mathbf{X}_0 = \left( X_0, ..., X_{-d+1} \right)^{\prime}$, that is, a distribution which results in a strictly stationary process $( X_t )_{ t \geq - d + 1 }$ , if and only if all (possible complex) roots of the equation

have modules strictly less than one. Then, the quadratic characteristic of the martingale difference array $\left( \xi_{nt}, \mathcal{F}_{nt} \right)$ satisfies

in probability as $n \to \infty$, since $\mathbf{X}_{t-1}$ is measurable w.r.t $\mathcal{F }_{n, t-1}$ and $\epsilon_t$ and $\mathcal{F }_{n, t-1}$ are independent. Moreover, the Lindeberg condition is satisfied, and the multivariate CLT implies that

Therefore, we obtain the CLT for the least squares estimator $\hat{\rho}_n$ as below

Consider also the $\left\lVert \ . \ \right\rVert_{ \Sigma }$ to be the norm on $\mathbb{R}^d$ pertaining to the inner product induced by $\Sigma$, that is, $\left\lVert x \right\rVert_{ \Sigma } = \left\lVert \Sigma^{1/ 2} x \right\rVert$ for all $x \in \mathbb{R}^d$, where $\left\lVert . \right\rVert$ is the Euclidean norm in $\mathbb{R}^d$.

Then, for any sequence $x_n \to \infty$ with $x_n = \mathcal{O}( n^{\delta - \epsilon} )$ for some $\epsilon > 0$ as $n \to \infty$, we have the moderate deviations limit theorem

Consider the martingale difference array such that $\tilde{\xi}_{nt} = \sigma^{-2} \Sigma^{-1 / 2} \xi_{nt}$. Using this martingale difference array we have that

Estimators for Gaussian Autoregressive Process

Consider the OU diffusion solution of the stochastic differential equation

where $W$ is a standard Brownian motion and $\theta$ is an unknown parameter in $\mathbb{R}$. For a continuous observation of $X$ over $[0,T]$, it is usual to consider an estimating function, the score function, that is, the derivative of the log-likelihood, which is given by

We denote by $\hat{ \theta }_T$ the maximum likelihood estimator (MLE), solution of $Y_T ( \theta ) = 0$, explicitly as

Thus, as in yu2009moderate, for estimating $\theta$, the following two estimators are widely used

Consider the least-squares estimator $\theta_n$ of the parameter $\theta$ which is defined as below (see, worms2001large)

It is well-known that $\left( \hat{\theta}_n \right)$ is strongly consistent, but the behaviour in distribution and the corresponding speeds are different according to the true parameter value of $\theta$, when $X_0 = 0$.

The Exogeneity Assumption in Time Series Models

The second important issue when considering econometric modelling is the aspect of exogeneity. In other words, the concept of exogeneity is essential to uniquely characterize the implications that certain variables are exogenous according to particular definitions.

The concept of structurally invariant conditional models characterizes the conditions which guarantee the appropriateness of "policy implications" or other control exercises, since any change in the distribution of the conditioning variables has no effect on the conditional model and thus on the conditional forecasts of the endogenous variables (see, ericsson1998exogeneity, white2014granger). It might be more helpful to understand what is not causality, thus we consider the definition of Granger Noncausality.

Autoregressive Processes with Exogenous Variables

Following the framework of damon2005estimation, we consider the following autoregressive Hilbertian with exogenous variables of order one model, denoted by ARHX(1) such that

where $X_{t,1} + ... + X_{t,p}$ are $p$ autoregressive of order one exogenous variables associated respectively with operators $\varphi_1,..., \varphi_p$ and strong white noises $u_{t,1},..., u_{t,p}$ such that

Furthermore, define with

In particular, when the sequence $\left\{ Y_t \right\}$ is an ARHX (1) model, then $\left\{ \mathcal{Y}_t \right\}$ is an $H^{p + 1}-$valued ARHX(1) process, as seen by the following autoregressive representation

We suppose that the noises $\varepsilon_t$ and $u_{t,1},..., u_{t,p}$ are independent which ensures that for all $i \in \left\{ 1,..., p \right\}$ the whole processes $( \varepsilon_t )$ and $\left( X_{t,i} \right)$ are independent since $X_{t,i}$ can be expressed as an infinite moving average of the $u_{t-q,1}$ where $q = 0,..., \infty$.

Consider for example the framework given by seo2019cointegrated. A cointegrated linear process in Bayes Hilbert space is isomorphic to a cointegrated linear process in a Hilbert space of centered square-integrable real functions. We illustrate the use of this isomorphism for modeling nonstationary time series of probability densities. In particular, a recent literature on functional data analysis deals with datasets whose observations take values in an infinite dimensional Banach or Hilbert space, typically a space of functions.

Temporal Dependence

Econometric models provide parsimonious representations of economic phenomena; thus the development of robust statistical methods is of paramount importance in improving our understanding of financial markets and economic decision making. In this section we study some aspects of estimation and inference in time series regression under temporal dependence.

Setting the time series framework

Consider the following model

where $\left( X_{t} \right)_{ t \in \mathbb{Z} }$ is the stochastic time series we aim to model, $\left( M_{t} \right)_{ t \in \mathbb{Z} }$ is a covariate process and $\left( \epsilon_{t} \right)_{ t \in \mathbb{Z} }$ a noise process. Here, we notice that $X_t = f_t \left( X_{t-1} \right)$ holds for the random function defined by $f_t (x) = F \left( x, M_{t-1}, \epsilon_t \right)$. Therefore, the stochastic time series model (ref) shows that the stochastic sequence $\left\{ X_t \right\}$ transforms a sequence of i.i.d random innovations $\epsilon_t$ into a time series sequence of dependent random variables (maps) as described by the random function $\left( f_t \right)_{t \in \mathbb{Z} }$. Moreover, the stochastic process $\left( M_{t} \right)_{ t \in \mathbb{Z} }$ allows to incorporate in the model exogenous covariates which typically exhibit temporal dependence. Under this general setting we can examine both linear and nonlinear time series models without restricting the estimation methodology employed.

Example: Nonlinear time series models

A nonlinear time series methodology which is commonly employed to capture conditional heteroscedasticity and dynamic variances in financial markets includes the ARCH and GARCH models. Specifically, these two nonlinear type of models have been the cornerstone of financial econometrics and statistical risk management fields. The ARCH model has been generalized to to GARCH by bollerslev1986generalized. A GARCH$(p,q)$ sequence $\left\{ X_t, - \infty < t < \infty \right\}$ has the following form

\color{black}

The literature has evolved considerably since the introduction of the GARCH$(p,q)$ model. The aforementioned stochastic linear and nonlinear time series models exhibit different asymptotic behaviour under the presence of exogenous covariates, especially considering these time series models as iterations of dependent random variables capturing linear and nonlinear dynamics.

Furthermore, considering these time series models as stochastically recursive sequences, when the independence assumption is removed under certain regulatory conditions, certain interesting asymptotic cases appear. For instance, we can study stationarity of asymptotic power GARCH processes which encompasses the case of exogenous regressors in nonlinear dynamics. Notice that for strictly exogenous regressors, the processes $\left( M_{t} \right)_{ t \in \mathbb{Z} }$ and $\left( \epsilon_{t} \right)_{ t \in \mathbb{Z} }$ are independent.

We are interested for example in the dynamics of predictive regression models described by (ref) with regressors that are not necessarily strictly exogenous, assuming that at any time $t$, the noise $\epsilon_t$ is independent from the past information $\sigma \left( \left( M_s, \epsilon_s \right) : s \leq t - 1 \right)$. Therefore, this independence assumption is weaker than the independence assumption between the two processes $\left( \epsilon_{t} \right)_{ t \in \mathbb{Z} }$ and $\left( M_{t} \right)_{ t \in \mathbb{Z} }$. The latter independence condition implies strict exogeneity, a terminology intially defined by Sims. Strict exogeneity is useful for deriving the conditional likelihood of the $X_t$'s conditionally on the $M_t$'s. However, strict exogeneity is a strong assumption. Moreover, chamberlain1982multivariate shown has shown that this assumption is equivalent to the non-Granger causality, that is, $M_t$ is independent of $\left( X_{s} \right)_{ s \leq t }$ conditionally on $\left( X_{s} \right)_{ s \leq t - 1 }$. This idea means that the covariate process $\left( M_{t} \right)_{ t \in \mathbb{Z} }$ evolves in a totally autonomous way. On the other hand, our exogeneity condition allows general covariates of the form $M_t = H \left( \eta_t, \eta_{t-1}, ...., \right)$ with $H$ a measurable function and a sequence $\left( \eta_t, \epsilon_t \right)_{t \in \mathbb{Z}}$ of i.i.d random vectors, $\epsilon_t$ being possibly correlated with $\eta_t$. Thus, the error $\epsilon_t$ can then still have an influence on future values of the covariates. For linear models, the two technical independence conditions discussed above between the noise and the covariate processes are often used as a distinction between weak and strict exogeneity (see, debaly2021iterations).

The family of ARCH and GARCH models are commonly used parametric models for capturing important stylized features of financial time series such as long memory, volatility clustering, fat tails in the distribution of stock returns etc. In particular, the GARCH(1,1) model or AR(1)-GARCH(1,1) model have repeatedly proved to provide robust parsimonious representations of market volatility (see berkes2003garch). A GARCH(1,1) regression process is given by the following equations

An AR(1)-GARCH(1,1) regression process is given by

The above non-linear models are valid under certain regulatory conditions which include $\omega , \alpha, \beta \geq 0$ and $\alpha + \beta < 1$ in order to avoid the existence of IGARCH effects. The recursively estimated sequences $(\sigma^2_t)$ for $t=0 ,..., + \infty$ are assumed to be non negative with probability 1 and common unconditional local mean given by $\sigma^2=\frac{\omega}{1 - \alpha - \beta}$. Statistical estimation can be done by imposing certain parametric assumptions. A commonly used methodology in the literature is the implementation of the QMLE.

The use of QML estimation can restore any inefficient model estimations (i.e., existence of negative coefficients or out of bounds persistence) and allows for higher estimation precision since the proposed methodology is less restrictive in the moment assumptions of the observed process (no restriction on normality assumption - which induces distribution-free inference).

and the QMLE for the set of model parameters $\theta=\{\omega, \alpha, \beta \}$ of GARCH(1,1) to be

which is equivalent to

Then asymptotic normality holds, which implies that,

Heavy Tail Time series

Max-stable processes

Tail behaviour of sample quantiles

Moreover, we have that

with the sum converging a.s for every $t \in T$. We have already proved that each $\left( X_j (t), t \in T \right)$ has a measurable version. Since the pointwise limit of measurable functions is measurable, we have constructed our measurable version of $\left( X (t), t \in T \right)$.

Extreme Value Theory for nonstationary time series

Consider the limit theory for extreme values of a class of nonstationary time series with the form

where $\left\{ Z_t = \sigma_t \eta_t , - \infty < t < \infty \right\}$ and $\left\{ \eta_t ; - \infty < t < \infty \right\}$ is a sequence of i.i.d random variables with regularly varying tail probabilities. Some convergence results for point processes based on one-sided moving averages $\left\{ \xi_t = \sum_{j=0}^{\infty} c_j Z_{t-j} \right\}$ are derived. Extreme properties of the nonstationary sequence $\left\{ Y_t \right\}$ are then obtained from the convergence results (see, niu1997extreme).

We assume that

and

where $0 \leq \pi_0 \leq 1$. Consider one-sided moving averages of the form

where $\left\{ c_j \right\}$ is a sequence of real constants with $c_0 = 1$

where

Furthermore, we assume that

A linear process can be written in the following form

is regularly varying with index $\alpha > 0$ if the i.i.d sequence $\left( \epsilon_{t} \right)$ is regularly varying with index $\alpha$. Under mildly conditions on $\left( \psi_j \right)$ we have that

For an i.i.d sequence $\left( A_t, B_t \right)_{t \in \mathbb{Z}}$, $A, B > 0$, the stochastic recurrence equation can be written as below

has a unique stationary solution given by the following expression

Limit Theory under Network Dependence

Random Graphs and Moderate Deviations Theory

Spectral Edge in Sparse Random Graphs: Tail Large Deviations

Understanding asymptotic properties of spectral statistics arising from random matrices has been the subject of intense study in recent years. Within this theme, a particularly important research direction is in deriving large deviation principles (rare event probabilities) for spectral functionals, such as the empirical spectral measure and the extreme eigenvalues. In the results above, we consider the large deviations of the edge eigenvalues of $\mathcal{G}_{n,p}$, in the regime $p << \sqrt{\text{log} n / \text{log} \text{log} n}/ n$. Notice for example that the largest eigenvalue $\lambda_1 \left( \mathcal{G}_{n,p} \right)$ is asymptotically $\sqrt{L}_p$, which is governed by the maximum degree of the graph, and corresponding eigenvector is localized; whereas above this threshold the typical value is asymptotically $np$, which is dictated by the total number of edges in the graph, and the corresponding eigenvector is completely delocalized.

Limit theorems for network dependent processes

Network dependence condition

We discuss in details the framework proposed by kojevnikov2021limit as this is useful in understanding the relevant limit theory. In particular, we provide a sufficient condition for the shape of the network that ensures that our limit theorems hold (LLN and CLT). A crucial aspect for the network based limit theorems to hold, the number of neighbours at distance $s$ should not grow too fast as $s$ increases. The precise condition for such neighbourhood shells depends on the dependence coefficients $\theta_{n,s}$, so that if $\theta_{n,s}$ decreases fast as $s$ increases, the requirement for the neighbourhood shells can be weakened. To introduce sufficient conditions, let

where $N_n^{\partial} (i;s)$. When $k = 1$, we simply write that $\delta_n^{ \partial } (s;1) = \delta^{ \partial }_n(s)$. This quantity measures the denseness of a network. Let us introduce a further notation. Define the following expression

where $N_n^{\partial}(i;s)$ is defined above, and we take $N_n(j;s-1) = 0$ if $s = 0$. We also define the following expression

The quantity $c_n (s, m; k)$ is easy to compute when a network is given, and it captures the network properties that are relevant for the limit theorems. It consists of two components: $\Delta_n (s,m; k \alpha )$ and $\delta_n^{\partial} \left( s; \alpha / (\alpha - 1) \right)$. These capture the denseness of the network through the average neighbourhood sizes and the average neighbourhood shell size. We summarize a sufficient condition for the network and the weak dependence coefficient as follows

Below we show that Condition ND (i) is sufficient for the LLN and CLT. Notice that $\Delta_n (s , m_n ; k)$ tends to decrease fast to zero as $s$ goes beyond a certain level, because the set $N_n ( j ; s-1)$ quickly becomes large.

Law of Large numbers

Let $\left\{ Y_{n,i} \right\}$ be conditionally $\psi-$dependent given $\left\{ \mathcal{C}_n \right\}$. Since a LLN can be applied element-by-element in the vector case, without loss of generality we can assume that $Y_{n,i} \in \mathbb{R}$ in this sections, that is, $v = 1$.

Define with,

Assume the following moment condition holds.

The next assumption puts a restriction on the denseness of the network and the rate of decay of dependence with the network distance.

For example, Assumption 4 above can fail, for example, if there is a node connected to almost every other node in the network as in the following example. Consider a network with the star topology, which has a central node or hub connected to every other node.

An unconditional version of the result, which replaces the conditional norm in Theorem 3.1 with the unconditional norm, can be established in a similar manner by replacing the conditional moment with the unconditional moment.

Next, we discuss LNNs for nonlinear functions of $\left\{ Y_{n,i} \right\}$. When $f \in \mathcal{L}_{ \nu, 1}$, a LLN for a nonlinear transformation $f \left( Y_{n,i} \right)$ follows immediately from the definition of the $\psi-$ dependence.

In that case, we have that

We have the following result.

However, in general nonlinear transformations of $\psi-$dependent processes are not necessarily $\psi-$dependent. In such cases, LLNs for nonlinear transformations can be established using the covariance inequalities for transformation functions. For example, suppose that for some nonlinear function $h(.)$ of a $\psi-$dependent process $\left\{ Y_{n,i} \right\}$, and that $\theta_{n,s}$ is bounded by a constant uniformly over $s \geq 1$ and $n \geq 1$.

In that case for some constants $C > 0$ and $p > 2$, the conditional covariance given $\mathcal{C}_n$ between $\big\{ h(Y_i) - \mathbf{E} \left[ h( Y_{n,i} ) | \mathcal{C}_n \right] \big\}$ and $\big\{ h( Y_{n,j} ) - \mathbf{E} \left[ h( Y_{n,j} ) | \mathcal{C}_n \right] \big\}$ is bounded by

Therefore, as $n \to \infty$,

provided that sup$_{n,i} \left\lVert h \left( Y_{n,i} \right)\right\rVert_{ \mathcal{C}_{n,p}} < \infty$ almost surely.

Central limit theorem

In this section, we study the CLT for a sum of random variables that are conditionally $\psi-$dependent. Define with

It is not hard to see that Condition ND is a sufficient condition for this assumption, when $\sigma_n \geq c \sqrt{n}$ with probability one, for some constant $c > 0$ that does not depend on $n$. The latter condition is satisfied if the "long-run variance", Var $\left( S_n | \mathcal{C}_n \right) / n$ is bounded away from $c^2 > 0$ for all $n \geq 1$. The theorem below establishes the CLT for the normalized sum $S_n / \sigma_n$.

The proof of the CLT uses Stein's Lemma and the CLT immediately gives a stable convergence of a normalized sum of random variables under appropriate conditions. Suppose that

where $v^2$ is a random variable that is $\mathcal{C}-$measurable and $\mathcal{C}$ is a sub $\sigma-$field of $\mathcal{C}_n$ for all $n \geq 1$. Then, it follows that $S_n / \sqrt{n}$ converges stably to a mixture normal random variable.

Network HAC

In this section, we develop network HAC estimation of the conditional variance of $S_n / \sqrt{n}$ given $\mathcal{C}_n$, where $S_n = \sum_{ i \in N_n } Y_{n,i}$. We assume that $\mathbf{E} \left[ Y_{n,i} | \mathcal{C}_n \right] = 0$ almost surely for all $i \in N_n$. Let

Then, the conditional variance of $S_n / \sqrt{n}$ given $\mathcal{C}_n$ is expressed as below

Similarly to the time-series case, the asymptotic consistency of an estimator $V_n$ requires a restriction on weights given to the estimated "autocovariance" terms $\Omega_n(.)$. Consider a kernel function $\omega: \mathbb{R} \to [-1,1]$ such that $\omega(0) = 1$, $\omega(z) = 0$ for $|z| > 1$, and $\omega(z) = \omega( -z)$ for all $z \in \mathbb{R}$. Let $b_n$ denote the bandwidth or the lag truncation parameter. Then, the kernel HAC estimator of $V_n$ is given by

where $\omega_n(s) = \omega( s / b_n )$, and

We have that $\bar{Y} = S_n / n$ is a consistent estimator of $\Lambda_n$ such that

We redefine the kernel HAC estimator as follows

where

Stable Limit Theorems under Network Dependence

We borrow some of the derivations presented in the framework of lee2019stable. A relevant question of concern is: Does imposing conditions on the degree of connectedness in the network allows to control the stability of the nodes?

Stable Convergence of an Empirical Process

Suppose that $\left\{ Y_i \right\}_{ i \in N_n }$ is a given triangular array of $\mathbb{R}-$valued random variables with is CND with respect to $\left( \nu_n, \mathcal{M} \right)$. Let $\mathcal{H}$ be a given class of real measurable functions on $\mathbb{R}$, having a measurable envelope H. Then, we consider the following empirical process:

where for each $h \in \mathcal{H}$, we have that

The empirical process $\nu_n$ takes a value in $\ell^{\infty} \left( \mathcal{H} \right)$, the collection of bounded functions on $\mathcal{H}$ which is endowed with the sup norm so that $\big( \ell^{\infty} \left( \mathcal{H} \right), \left\lVert . \right\rVert_{\infty} \big)$ forms the metric space $\left( \mathbb{D}, d \right)$ with the sup norm $\left\lVerth\right\rVert_{\infty} \equiv \text{sup}_{y \in \mathbb{R}} \left| h(y) \right|$. lee2019stable explore conditions for the class $\mathcal{H}$ and the joint distribution of the triangular array $\left\{ Y_i \right\}_{ i \in N_n }$ which delivers the stable convergence of the empirical process. Stable convergence in completer separable metric spaces can be defined as a weak convergence of Markov kernels. However, this definition does not extend to the case of empirical processes taking values in $\mathbb{D}$ that is endowed with the sup norm, due to the non-measurability. Weak convergence of an empirical process to a Gaussian process is often established in three steps. First, they show that the class of functions is totally bounded with respect to a certain pseudo-metric $\rho$. Second, they show that each finite dimensional projection of the empirical process converges in distribution to a multivariate normal random vector. Third, they establish the asymptotic $\rho-$equicontinuity of the empirical process.

Maximal Inequality

This subsection presents a maximal inequality in terms of bracketing entropy bounds. The maximal inequality is useful primarily for establishing asymptotic $\rho-$equicontinuity of the empirical process. We begin with a tail bound for a sum of CND random variables. The following exponential tail bound is crucial for our maximal inequality.

Furthermore, if Condition A holds and the condition $\mathbf{E} \left[ X_i | \mathcal{M}_{\nu_n(i)} \right] = 0$ is replaced by $\mathbf{E} \left[ X_i | \mathcal{G} \right] = 0$ and the $\sigma-$fields $\mathcal{M}_{\nu)n(i)}$ in $\sigma_i'$s are replace by $\mathcal{G}$, then the following hods for any $\eta > 0$,

Stable Central Limit Theorem

Moderate Deviations: Eigenvalues of Covariance Matrices

Although due to the convenient structure of the covariance matrix, the estimation procedure of the entries of the matrix is many times overlooked therefore the spectral analysis of high dimensional matrices mainly focused on the asymptotic behaviour of the eigenvalues without taking into consideration how the entries are estimated. In our framework we give emphasis on the estimation of the risk matrix with the use of quantile predictive regression models which clearly accounts for more features regarding the stationarity and persistence properties of the time series. In practise, the estimation of the covariance matrix with the conventional methodologies which consider the use of a high dimensional vector does not account for much features of the time series under examination (see, katsouris2021optimal, katsouris2023statistical).

According to zhang2020estimation, in recent years, large sample properties for high-dimensional sample covariance matrices, including their eigenvalues and eigenvectors, have been proved to be useful. In fact, random matrix theory (RMT) provides a plethora of useful methods for estimation and testing in high dimensional environments. Currently, there are two steams of literature about asymptotic theory of the largest eigenvalues of high dimensional random matrices. The first stream of literature is concerned with the Tracey-Widom law. In particular, it is well known that the limiting distribution of the largest eigenvalues of high-dimensional random matrices, such as Wigner matrices, follow the Tracy-Widom law which holds for Gaussian Wigner ensembles. Moreover, the second stream of literature is interested in the asymptotic behaviour of spiked eigenvalues.

Most of the existing studies operate under the assumption that observations of high dimensional data are independent. However, applications in finance and economics with high dimensional settings involve data which can be temporally dependent or even nonstationary. For instance, the first framework that examines the asymptotic behaviour of the largest eigenvalues of sample covariance matrices generated from high-dimensional nonstationary time series is proposed by zhang2018clt. Investigating the asymptotic behaviour of the largest eigenvalues of random matrix models is important in understanding the stability behaviour of these systems. See Remark 4 (page 34) of grama2006asymptotic who mention: "The key point in the proof when considering the eigenvalues of a covariance or state matrix, is the fact that the eigenvalues of a matrix depend continuously on its elements. In practise, this means that the maximal eigenvalue of a matrix is a predictable discrete time process in $t$ and therefore the random times $\tau_m$ are predictable stopping times". \color{black}

Limit of the smallest eigenvalue of a Large Sample Covariance Matrix

The related framework is presented by bai1999exact.

\paragraph{Proof of Lemma 3 of the paper}

The following inequality holds,

Limit Theory for the largest eigenvalues of sample covariance matrices with heavy-tails

Theory that consistently estimates the spectrum of a large dimensional covariance matrix using RMT. Thus, statistical considerations will be our motivation for a random matrix model with heavy-tailed and dependent entries. Notice that results on the global behaviour of the eigenvalues of $XX^{\top}$ mostly concern the spectral distribution, that is, the random probability measure of its eigenvalues $p^{-1} \sum_{i=1}^p \epsilon_{ n^{-1} \lambda_{(i)} }$, where $\epsilon$ denotes the Dirac measure. The spectral distribution converges, as $n, p \to \infty$ with $p / n \to \gamma \in (0,1]$, to a deterministic measure with density function

This is the so-called Marcenko-Pastur law. Therefore, one obtains a different result if $X X^{\top}$ is perturbed via an affine transformation. Although the eigenvalues of $XX^{\top}$ offer various interesting local properties to be examined, we will only focus on the joint asymptotic behaviour of the $k$ largest eigenvalues $\left( \lambda_{(1)}, ..., \lambda_{(k)} \right), k \in \mathbb{N}$. This is motivated from a statistical point of view since the variances of the first $k$ principal components are given by the $k$ largest eigenvalues of the covariance matrix. Therefore, a Marcenko-Pastur type result holds for the asymptotics of the spectrum of heavy-tailed random matrices.

Limit Spectral Distribution for symmetric random matrices with correlated entries

Let $\left( X_{k, \ell} \right)_{ ( k, \ell ) \in \mathbb{Z}^2 }$ be an array of real-valued random variables, and consider its associated systemic random matrix $\mathbf{X}_n$ of order $n$ defined by

We then define

The aim of this section is to study the limiting spectral empirical distribution function of the symmetric matrix $\mathbf{X}_n$ when the process $\left( X_{k, \ell} \right)_{ ( k, \ell ) \in \mathbb{Z}^2 }$ has the following dependence structure: for any $(k , \ell) \in \mathbb{Z}^2$,

where $\left( \xi_{ i, j } \right)_{ (i,j) \in \mathbb{Z}^2}$ is an array of i.i.d real-valued random variables given on a common probability space $\left( \Omega, \mathcal{F}, \mathbb{P} \right)$ and $g$ is a measurable function $\mathbb{R}^{ \mathbb{Z}^2 } \mapsto \mathbb{R}$ such that $\mathbb{E} \left( X_{0,0} \right) = 0$ and $\left\lVert X_{0,0} \right\rVert_2 < \infty$.

The Theorem below shows a universality scheme for the random matrix $\mathbb{X}_n$ when the entries of the symmetric matrix $\sqrt{n} \mathbb{X}_n$ have the above dependence structure. Notice that the particular result does not require rate of convergence to zero of the correlation between the entries of the risk matrix.

The above theorem is important since it shows that the study of the limiting spectral distribution function of a symmetric matrix whose entries are functions of $\textit{i.i.d}$ random variables can be reduced to studying the same problem as for a Gaussian matrix with the same covariance structure.

\paragraph{Proof of Theorem (ref)}

Moderate Deviations for Extreme Eigenvalues

According to jiang2021moderate the Marcehnko-Pastur Law implies that the empirical spectral distribution of $W$ converges to a deterministic distribution with support $\left[ \left( 1 - \sqrt{ \beta } \right)_{+}^2 , \left( 1 + \sqrt{ \beta } \right)^2 \right]$, where $x_{+} = \text{max} \left\{ 0, x \right\}$.

As a corollary we can obtain the following moderate deviation for the condition number $\lambda_{\text{max} } / \lambda_{\text{min} }$, which is important quantity that can indicate whether the matrix is ill-conditioned.

Moderate Deviations for $\lambda_{ \text{max} }$, $\lambda_{ \text{min} }$ with fixed $k$

Suppose there exist two sequences $b_n \left( b_n \to \infty \right)$, $r_n$ so that a family of random variables $\left\{ Z_n , n \geq 0 \right\}$ with values in topological vector space $\mathcal{X}$ (equipped with $\sigma-$field $\mathcal{B}$) satisfies a fluctuation theorem, say $b_n ( Z_n - r_n )$ converges to some nontrivial distribution. We say $Z_n$ satisfies the moderate deviation principle with speed $\lambda_n \to \infty$ and with good rate function $I(.)$ if the level sets $\left\{ I \leq \ell \right\}$ are compact for all $\ell > 0$ and for any closed set $F$ and open set $G \in \mathcal{B}$,

Here, $\tilde{b}_n$ is a sequence satisfying $\tilde{b}_n \to \infty$ and $\tilde{b}_n / b_n \to 0$. The form of rate function $I(.)$ is closely related to the limit distribution of $b_n ( Z_n - r_n )$. Therefore, if the rate function is continuous we have that

Discussion on High Dimensional Results

The literature has indeed documented the inconsistency of sample covariance matrices for especially high-dimensional settings. In order to be able to investigate the features of our novel tail risk matrix, it is necessary to consider the distribution theory of the entries of the matrix. Overall to assess the quality of a matrix estimate, we can use the norm operator such as the spectral radius. Below we present some useful results (e.g. see, chen2013covariance).

Suppose we have $n$ temporally observed $p-$dimensional vectors $\left( \mathbf{z}_i \right)_{ i = 1}^n$ having mean zero and covariance matrix $\Sigma_i = \mathbb{E} \left( \mathbf{z}_i, \mathbf{z}_i \right)$ whose dimension is $p \times p$. Our goal is to estimate the covariance matrices $\Sigma_i$ and their inverses $\Omega_i = \Sigma_i^{-1}$ based on the data matrix $Z_{ p \times n} = \left( \mathbf{z}_1,..., \mathbf{z}_n \right)$. In the classical situation where $p$ is fixed, $n \to \infty$ and $\mathbf{z}_i$ are mean zero independent and identically distributed $\textit{i.i.d}$ random vectors.

In particular, it is well known that the sample covariance matrix $\hat{ \boldsymbol{\Sigma} }_n = \frac{1}{n} \sum_{ i = 1 }^n \boldsymbol{X}_i \boldsymbol{X}_i^{ \top }$, is a consistent and well behaved estimator of $\boldsymbol{\Sigma}$, and $\hat{ \boldsymbol{\Omega} }_n = \hat{\boldsymbol{\Sigma}}_n^{-1}$ is a natural and good estimator of $\boldsymbol{\Omega}$. However, when the dimensionality $p$ grows with $n$, then $\hat{\boldsymbol{\Sigma}}_n$ is no longer a consistent estimate of $\boldsymbol{\Sigma}$ in the sense that its eigenvalues do not converge to those of $\boldsymbol{\Sigma}$. Moreover, it is clear that $\hat{ \boldsymbol{\Omega} }_n$ is not defined when $\hat{\boldsymbol{\Sigma}}_n$ is not invertible in the high-dimensional case with $p >> n$. Define with

the eigen-decomposition, where $\mathbf{Q}$ is an orthonormal matrix and $\hat{ \boldsymbol{\Lambda} }$ is a diagonal matrix. For $v > 0$,

where $0 < v \leq \sqrt{p} \bar{\omega}$ and $\omega^2$ is the rate of convergence. Let $\mu_1,...,\mu_p$ be the diagonal elements of $\mathbf{Q}^{\top} \boldsymbol{\Sigma} \mathbf{Q}$. Then, we have by Theorem 2.1 that $\sum_{j=1}^p \left( \hat{\lambda}_j - \mu_j \right)^2 \leq p^2 \bar{\omega}^2$, and consequently

If $\hat{\lambda}_j \leq 0$, since $\mu_i \geq 0$, we have that $\left| \hat{\lambda}_j \right| \leq \left| \hat{\lambda}_j - \mu_i \right|$. Then,

The eigenvalues of $\tilde{S}_v$ are bounded below by $v$, and thus it is positive definite. Suppose that

Thus, the same positive-definization procedure also applies to the spectral norm. We obtain that the differences between the eigenvalues and eigenvectors of $\Sigma_m$ and $\hat{\Sigma}_m$ can be bounded by the following

Moreover, we have that

with

where we denote with $S_{r,m} = \sum_{ s \neq r } \frac{1}{ \lambda_s - \lambda_r } e_{s,m} e_{s,m}^{\top}$. Assumption 1 implies that $\mathbb{E} \left( \hat{ \beta }_r \right) = 0$ and Var$\left( \hat{ \beta }_r \right) = \frac{ \lambda_r }{ r }$ and with $\delta_{ii} = 1$ as well as $\delta_{ij} = 0$ for $i \neq j$, we obtain that

for all $m$.

Furthermore, since

Therefore after applying the central limit theorem we obtain that

Obviously the event $\hat{\lambda}_{r-1} > \hat{\lambda}_{r} > \hat{\lambda}_{r+1}$ occurs with probability 1.

\color{black}

Nonstationary Time Series Regressions

Testing for unit root in time series regression

Consider the sequence of partial sums $\left\{ S_t \right\}$, then we obtain that

for $\frac{( j - 1)}{T} \leq r \leq \frac{j}{T}$ for $j \in \left\{1,..., T \right\}$, where $[Tr]$ denotes the integral part of Tr. Then, $X_T(r)$ lies in $\mathcal{D} = \mathcal{D}[0,1]$, the space of real valued functions on the interval $[0,1]$ that are right continuous and have finite left limits. Under very general conditions the random element $X_T(r)$ obeys a central limit theory on the function space $\mathcal{D}$. In particular, we have that, as $T \to \infty$,

which implies a weak convergence on the associated probability measure. In particular, in this case, the probability measure of $X_T(r)$ converges weakly to the probability measure of the standard Brownian motion $W(r)$. We focus on the asymptotic behaviour of the sample moments of the process $\left\{ S_t \right\}$ and the innovations $\left\{ u_t \right\}$. We express the limit distributions as functions of standard Brownian motion $W(r)$.

All integrals are understood to be taken over the interval $[0,1]$, while integrals such as $\int W$, $\int W^2$, $\int r W$ are understood to be taken with respect to Lebesgue measure such that we write $W_1 = W(1)$. Then, as $T \to \infty$, we obtain that

Power Functions for Unit Root Tests

Consider the power functions for unit root tests by considering the sequence of local alternatives as

where $c=0$ reduces to the null hypothesis, $c > 0$ gives local explosive alternatives and $c < 0$ corresponds to local stationary alternatives. Moreover, the asymptotic theory for the sample moments of near-integrated time series converge weakly to corresponding functionals of a diffusion process rather than standard Brownian motion. Specifically, we have that, as $T \to \infty$,

where

is the OU process generated in continuous time by the stochastic difference equation

A particular feature of the local to unity model is that the localizing parameter is identifiable but it is not consistently estimable.

Consider the following expression which gives the signal to noise ratio

and so the signal from $x_t$ is too weak relative to the error variation to produce a consistent estimator of the localizing coefficient $c$. Although methods have been developed to utilize the way in which the limit distribution depends on the localizing coefficient, the failure of a consistent estimation has been a challenging task when conducting inference in these models. The dependence of the limit distribution on $c$ also affects resampling procedures such as the bootstrap, which are known to be inconsistent in models of this type because exactly of this dependence of functionals to the nuisance parameter of persistence.

Time Series Regression with a Unit Root

The consistent estimates $s_u^2$ and $s^2_T$ are used to develop new tests for unit roots that apply under very general conditions. We define the statistics

where $Z_{ \alpha }$ is a transformation of the standardized estimator $T \left( \hat{\alpha} - 1 \right)$ and $Z_t$ is a transformation of the regression $t$ statistic. Then, the limiting distribution of $Z_{ \alpha }$ and $Z_t$ are given by

Regression Asymptotics using Martingale Convergence Methods

Considers weak convergence arguments to stochastic integral approximations. In particular for the case of a martingale-difference sequence $\left( \epsilon_t^2 | \mathcal{F}_{t-1} \right) = \sigma^2_{\epsilon}$ for all $t$ and $\mathsf{sup}_{ t \in \mathbb{Z} } \mathbb{E} \left( \epsilon_t^p \ \mathcal{F}_{t-1} \right) = \sigma^2_{\epsilon} < \infty$ almost surely for some $p > 2$. Then by Donsker's theorem for the partial sum process we have that

where $\big\{ W(s), s \geq 0 \big\}$ denotes the standard Brownian motion, which implies that the bilinear form:

Notice that the particular approach also has drawbacks. One of the drawbacks is that the approach is problem specific in certain ways. For instance it cannot be directly used in the case of statistics such as $\sum_{t=1}^n y_{t-1} u_t$, where $y_t = \alpha_n y_{t-1} + u_t$, for $t = 1,...,n$ and $\alpha_n \to 1$ as $n \to \infty$, that are central to the study of local deviations from a unit root in time series regression. Strong approximations to partial sums of independent random variables, along with an application of PS device, allow one to obtain invariance principles under independence and stationarity assumptions with explicit rates of convergence. For instance, using the Hungarian construction if $( \epsilon_t )_{t \in \mathbb{Z} }$ satisfy the conditions then

Notice that some existing results available in the literature on convergence to stochastic integrals can be applied to obtain convergence results. Denote with $ N_{n,r} = \frac{ \epsilon_0 }{ \sqrt{n} } + \sum_{t=1}^{ \floor{nr} } \frac{ \epsilon_t }{ \sqrt{n} }$. Then, the following stochastic integral representation holds for the statistic on the right-hand size

Since $N_{n,r} \to W(r)$, then the latter result implies that the convergence result to the stochastic integral holds, provided that the sequence of processes $\left\{ N_{n,t} \right\}$ satisfies the uniform tightness condition.

Unification of the Limit theory of Autoregression

Consider the following first order autoregression model (e.g., see, giraitis2006uniform)

with martingale-difference errors $\epsilon_t$. We treat the stationary $| \alpha | < 1$, unit root $\alpha = 1$, local to unity, and explosive cases together in what follows and show how the limit theory for all these cases may be formulated in a unified manner within the martingale convergence framework. Consider the stationary and unit root cases. For $r \in (0,1]$, define the recursive least squares estimator $\hat{\alpha}_r = \sum_{t=1}^{ \floor{nr} } y_{t-1} y_t / \sum_{t=1}^{ \floor{nr} } y_{t-1}^2$,

where $X_n(r)$ is the martingale given by

Then, by Theorem 4.1 it holds that

Then, it can be proved that

which unifies the limit theory for the stationary and unit root autoregression.

Defining the error variance estimator $s_r^2 = \frac{1}{ \floor{nr} } \sum_{ t=1 }^{ \floor{nr} } \big( y_t - \hat{\alpha}_r y_{t-1} \big)^2$ and noting that $s_r^2 \overset{ p }{ \to } \sigma^2_{\epsilon}$ for $r > 0$, we have the corresponding limit theory for the recursive t-statistic such that the following limit result holds

The theory also extends to cases where $\alpha$ lies in the neighborhood of unity such that $\alpha = \left( 1 + \frac{c}{n}\right)$. Then, it follows that

where $J_c (v) = \displaystyle \int_0^v e^{ c(v-s)} dW(s)$ is a linear diffusion. Then, we obtain that

Furthermore, we can investigate the effect on the asymptotic theory of the proposed test statistics when there are moderate deviations from unity of the form $\alpha = \left( 1 + \displaystyle \frac{c}{ n^{b} } \right)$ for some $b \in (0,1)$ and $c < 0$. In particular, under the assumption that the model covariates are generated as near unit root processes with $b \in (0,1)$ and $c < 0$ then the regressors are considered to be mildly integrated and the functional takes the following form

Time Series Regression

One of the main differences with independent data where the OLS estimator satisfies certain conditions is that, time-series data commonly violate one of the classical assumptions. In particular, the independence assumption holds that the error terms corresponding to different point in time are not correlated. However, when the error terms are serially correlated (autocorrelated), the OLS method produces biased estimates of the standard errors of the regression coefficients. The persistence of a random shock is what distinguishes an AR error process from an MA error process. An AR(p) error process can be written as an infinite sum of past random shocks. Each shock persists indefenitely, although its importance diminishes over time. For example, in an AR(1) process, the effect persists indefenitely, although it decreases over time and is effectively zero after some relatively small number of time periods ( choudhury1999understanding).

Unit Roots versus Deterministic Trends

Consider the following model

Then, assume that the true data generating process is given by

The Engle-Granger Approach

Consider a vector time series process $X_t \in \mathbb{R}^k$ such that each component is $X_{it}$ is integrated or order one, such that each of the series $X_{it}$ contains a unit root, but $\Delta X_{it}$ is a zero mean stationary process. Thus, using the Wold decomposition theorem we can write:

where $L$ is the lag operator and $u_t$ is a $k-$variate white-noise process.

Working in the context of a bivariate system with at most one cointegrated vector, engle1987co propose estimating the cointegrated vector $\alpha = \left( 1, \alpha_2 \right)$ by regressing the first component $X_{1t}$ of $X_t$ on the second component $X_{2t}$, using OLS, and then testing whether the OLS residuals $\hat{Z}_t$ have a unit root using the augmented Dickey-Fuller test. Therefore, if the test rejects the unit root hypothesis, and thus accepts the cointegration hypothesis, one may substitute for $\alpha^{\prime} X_{t-1}$ the OLS residual $\hat{Z}_{t-1}$ and estimate the parameter matrices in $A^{*} (L)$ by OLS, assuming that $A^{*}(L)$ is a finite-order lag polynomial. Notice that the above approach is only applicable if there is at most one cointegrated vector. Systems with dimension greater than two, however, may have multiple cointegrated vectors.

Testing the hypothesis of no serial correlation

In some cases we are interested to test the hypothesis that the series $Y_t$ is serially uncorrelated against the alternative that is serially correlated. The most appropriate test of this hypothesis is based on the OLS of an AR$(p)$ model. Consider the following model

with $e_t$ a MDS. Thus, in this model the series $Y_t$ is serially uncorrelated if the slope coefficients are all zero. For example, when the cointegrating vector is known and the covariates of the cointegration vector are strictly exogenous, the optimal test for the null of a unit root in the known cointegrating vector collapes to the usual unit root test for a single variable. From a testing perspective, an unknown cointegrating vector is a nuisance parameter that is identified only under the alternative hypothesis, so optimal tests against the alternative of cointegration take the form of a weighted average test where the weights are over the possible cointegrating vectors.

A Cointegrated VAR

A VAR has several equivalent representations that are valuable for understanding the interactions between exogeneity, cointegration and economic policy analysis. To start, the levels form of the $s-$th order Gaussian VAR for $x$ is

where $K$ is an $N \times N_0$ matrix of coefficients of the $N_0$ deterministic variables $q_t$.

Suppose that we have a vector $Y_t = \left[ y_{1t}, y_{2t},..., y_{nt} \right]^{\prime}$ that does not satisfy the conditions for stationarity. One way to achieve stationarity might be to model $\Delta y_t$, rather than $y_t$ itself. However, differencing can discard important information about the equilibrium relationships between the variables. This is because another way to achieve stationarity can be through linear combinations of the levels of the variables. Thus, if such linear combinations exist then we have cointegration and the variables are said to be cointegrated.

Cointegration has some important implications:

Identification, Estimation and Forecasting with Time Series Regressions:

Cointegrating Regression

Distribution of the OLS Estimate for a Special Case

Following, hamilton2020time let $y_{1t}$ be a scalar and $\boldsymbol{y}_{2t}$ be a $( p \times 1 )$ vector satisfying

In particular, if $y_{1t}$ and $\boldsymbol{y}_{2t}$ are both $I(1)$ but $u_{1t}$ and $\boldsymbol{u}_{2t}$ are $I(0)$, then $d \equiv ( p + 1 )$, the $d-$dimensional vector $\left( y_{1t}, \boldsymbol{y}_{2,t-1} \right)^{\prime}$ is cointegrated with the cointegrated relation given by expression (ref). Therefore, by considering the case that the innovation sequence of the system follows a Gaussian distribution, we have that the following holds

Notice that although we consider a Gaussian cointegrated system we assume that the explanatory variables $\boldsymbol{y}_{2,t}$ in (ref) are independent of the error term $u_{1t}$ for all $t$. Therefore, conditional on $\big( \boldsymbol{y}_{2,1}, \boldsymbol{y}_{2,2},..., \boldsymbol{y}_{2,T} \big)$, the OLS estimates have a Gaussian distribution given by

Therefore, the Wald form of the OLS $F-$test of the null hypothesis is given by

where

In other words, conditional on the vector $\big( \boldsymbol{y}_{2,1}, \boldsymbol{y}_{2,2},..., \boldsymbol{y}_{2,T} \big)$ it follows that the above expression when we replace with the population variance, has a $\chi^2_m$ distribution.

Thus, conditional on the data $\big( \boldsymbol{y}_{2,1}, \boldsymbol{y}_{2,2},..., \boldsymbol{y}_{2,T} \big)$ the OLS F-test could be viewed as the ratio of a $\chi^2_m$ variable to the independent $\chi^2_{ T - p }$ variable $( T - p ) \sigma_T^2 / \sigma_1^2$ with numerator and denominator each divided by its degree of freedom. Therefore, the OLS F test has an exact $F ( m, T - p )$ distribution. Hence, despite the $I(1)$ regressors and complications of cointegration, the correct approach is to is to estimate the model by OLS and use standard $t$ or $F$ statistics to test any hypotheses about the cointegrating vector.

Estimation of cointegrating vectors

We discuss estimation of $\beta$ with the following representation

from the observable vector $z_t = \left( x_t^{\prime}, y_t \right)^{\prime}$ introduced in the previous section. The setting is or cointegration we assume that holds requires that the unobservable process $e_t$ is $I(d_e)$, $d_e < d_y$, while we assume that $\beta$ is identified for $t=1,...,n$

For a generic column vector or scalar sequence $a_t$, $t = 1,...,n$ define the discrete Fourier transform below

With also a column vector or scalar sequence $b_t, t=1,...,n$ possibly identical to $a_t$, define the cross periodogram as below

Now denote by $\lambda_j = 2 \pi j / n$, for integer $j$, the Fourier frequencies, and define the averaged cross-periodogram. The case $m = \floor{ n / 2}$, where $\floor{.}$ denotes the integer part, is of particular interest, as we deduce that

the mean-corrected sample covariance, with $\bar{a} = n^{-1} \sum_{t=1}^n a_t$. We observe that $\widehat{F}_{ab} (m)$ represent the contributions from frequencies $[1, \lambda_m ]$ to the sample covariances above. We estimate $\beta$ by the frequency domain least squares (FDLS) statistic given by

Cointegrating Vector and Equilibrium Points

A typical cointegrating system has the following form

and require the coefficient matrix $A$ to have stable roots, then the new system is a conventional SEM. Furthermore, it is a triangular system since it is written in a reduced form. In particular, when $\Sigma$ is not block-diagonal, the MLE is obtained by using OLS on the augmented regression equation such that

Consider $\widehat{\boldsymbol{\beta}}$ to be the single equation OLS estimator of $\boldsymbol{\beta}$. Then, it holds that

Then, consider the fully modified OLS estimator given as below

whose asymptotics are as below

Then, the further modification for the endogeneity of $y_{2t}$ is required to remove the correlation between the Brownian motion $S_1$ and $S_2$.

Therefore, this is achieved by constructing

Then, the fully modified OLS estimator employs both the serial correlation and endogeneity corrections and is given by the following expression

Therefore, with these corrections the new estimator $\boldsymbol{\beta}^{+}$ has the same asymptotic behaviour as the full system MLE. Observed that it is a two-step estimator, and relies on the preliminary construction of $y_{1t}^{+}$ and $\widehat{\delta}^{+}$. Then, fully modified test statistics are based on $\boldsymbol{\beta}^{+}$ which can be constructed in the usual way. Therefore, we define the t-ratios as below

Unit Roots, Cointegration and Structural Breaks

Testing for Structural Change in Cointegrated Systems

Following the framework of seo1998tests, and by the invariance principle proposed by phillips1986multiple, it holds that

Therefore, we need to show that

For instance, if $\left\{ \Phi(L) u_t \right\}$ is uniformly square integrable, we can then apply the following result.

Therefore, we need to show that

where $\Phi_1 (L) = \left( \Phi(L) - \Phi(1) \right) / (1 - L)$.

Moreover, a key ingredient for deriving the limiting distribution of the Wald-type and LR-type statistics is consider the weak convergence of partial-sum functionals. We define weak convergence on the space $C\left( [0,1] \right)$ with respect to the uniform metric. Thus, we define weak convergence of the projected sequence

Related asymptotic moments for the particular framework are as below:

Therefore, under the null hypothesis of no structural change, the model can be estimated by existing methods. Then, it holds that

Notice that the expression $\boldsymbol{B}(s) = BM ( \boldsymbol{\Omega} )$ represents a Brownian motion with long-run variance $\boldsymbol{\Omega}$. Due to uniformly integrability, the following limit result holds

By Lemma 9 of seo1998tests, it holds that

uniformly in $\uptau \in \mathcal{T}$.

Thus, to derive the limit distribution under the null hypothesis, $ \frac{1}{n} \sum_{t=1}^{ \floor{\lambda n } } z_t z_t^{\prime} \lambda Q$, uniformly in $\lambda \in [0,1]$, for Q some positive definite matrix and for all $t$ and $n^{- 1 / 2} \sum_{t= 1}^{ \floor{ \lambda n } } z_t u_t \Rightarrow s Q^{- 1 / 2} W_q (s)$, where $W_q(s)$ is a $q-$vector of independent Wiener processes. This imposes a form of moment homogeneity, that is, a homogeneous distribution throughout the sample. The particular assumption holds only under the null of no structural break and is a necessary condition for deriving the limiting distributions of the test statistics.

Trend Stationarity and Structural Break Testing

Consider the following model for the scalar random variable $x_t$ such that

The shocks, $e_t$, are assumed to follow a zero mean fractionally integrated process of order $d$, denoted by $e_t \in I(d)$. Furthermore, the deterministic trend break term, $DT_t \left( \uptau^* \right)$, is defined for a generic $\uptau$ such that $DT_t \left( \uptau \right) := \left( t - \floor{\uptau T} \right) \boldsymbol{1} \left\{ t \geq \floor{\uptau T} \right\}$. Therefore, where a true trend break occurs, which implies that we are under the alternative hypothesis and so $\beta_3 \neq 0$, we assume that the true break fraction is such that $\uptau^* \in [ \uptau_U, \uptau_L ] =: \Lambda \subset [0,1]$, where the quantiles $\uptau_L$ and $\uptau_U$ are trimming parameters. The formulation of the model is obtained using the vectors such that $z_t(\tau):= \left( 1, t, DT_t(\uptau) \right)^{\prime}$ and $\beta := \left( \beta_1, \beta_2, \beta_3 \right)^{\prime}$. Therefore, we can obtain the OLS estimate of $\beta$ as below

Moreover, we define the corresponding de-trended residuals as below

where the particular sequence is obtained for $t = 1,...,T$.

Therefore, under the null hypothesis we have that the estimated $\eta_t$ are obtained by taking the corresponding fractional differences of these OLS de-trended residuals. In other words, if the true break fraction $\uptau^*$, was known then one would simply evaluate the Likelihood ratio for break testing, denoted with $LM (\uptau)$ at the location $\uptau = \uptau^*$, and thus the resulting test statistic is denoted by $LM (\uptau^* )$. However, most of the empirical applications in the literature do not have a priori information regarding the exact location of the break-point location within the full sample. Our focus, is on the case where $\uptau^{*}$ is unknown and so, our proposed test will be based on evaluating $LM (\uptau)$ at $\hat{\uptau}$, which is obtained as the minimum RSS estimate

whose exact form is determined according the exact value of $d_0$ being tested under the null hypothesis.

Specifically if $d_0$ lies in the region $( -0.5, 0.5 )$ then we estimate $\uptau^{*}$ using the levels of the data and testing the null hypothesis that the long memory parameter in the levels data is $d_0$, whereas if $d_0$ lies in the range $(0.5, 1.5)$ we instead estimate $\uptau^{*}$ using the first differences of the data and test the null hypothesis that the long memory parameter in the first difference data is $( d_0 - 1 )$ (see, arai2007testing).

Therefore, in order to be able to determine the limiting distribution of the test statistic, we will need to determine the large sample behaviour of LM evaluated at the estimated break-point location, denoted by $LM ( \hat{\uptau} )$, which is done by comparing it to the infeasible LM statistic. In other words, the large sample behaviour of $LM ( \hat{\uptau} )$ clearly depends on the large sample properties of the estimates $\hat{\uptau}$ and $\widehat{\beta} ( \widehat{\uptau} )$ evaluated at $\uptau = \widehat{\uptau}$. Now assume that the underline data generating process has a break under the alternative hypothesis, which implies that we can define the Pitman drift such that $H_c: \theta := \theta_T = c / \sqrt{T}$.

Then, for $d_0 \in \left( -0.5, 0.5 \right)$, define a generic element $\alpha$ and the diagonal matrix such that

Under the alternative hypothesis $H_c$, by setting $r = 0$ and $\alpha = \delta_0$, the FCLT holds such that

Cointegration Testing under Structural Breaks

Further to the two previous examples, in this section we consider cointegration testing under structural breaks (see, campos1996cointegration, johansen2000cointegration and kasparis2012dynamic). In order to test for structural breaks in cointegration models a residual-based test can be derived from single equation models. However, the system equation approach appears limited when a structural break in cointegrating vectors is incorporated in the conventional specification form due to the presence of parameter instability (see, saikkonen2006break). The literature proposes two different estimation and testing approaches, the first approach consists of the test for the null of cointegration while the second approach considers the test for the null of no cointegration (see, wagner2022residual). Specifically, the authors consider the Lagrange Multiplier (LM) test and derive its limiting distribution which is free of nuisance parameter dependencies except for the number of $I(1)$ regressors (degrees of freedom) and the location of the structural break.

To derive the asymptotic theory of the test statistic we derive an invariance principle and a corresponding functional central limit theorem to the partial-sum process $n^{ -1/2 } \sum_{t=1}^{ \floor{nr} } v_t$, which weakly converges to a two-dimensional standard Brownian motion $B(r) = \left( B_1(r), B_2(r) \right)^{\prime}$. In particular, the test proposed by shin1994residual is constructed based on the regression residual of $y_{1t}$ on 1 and $y_{2t}$, denoted by $\hat{v}_{1t}$.

Thus, the test statistic is given by

We show that $V_n$ diverges to infinity but $V_n / n$ weakly converges to a distribution that is positive almost surely. Denote with $x_t = \left( 1, y_{2t} \right)^{\prime}$ and $e_t = v_{1t} + \mu_2 d_{2t}$ such that $y_{1t} = \left( \mu_1, \beta \right) x_t + e_t$. The residual $\hat{v}_{1t}$ is

The following weak convergence results hold

Furthermore, it holds that

where the normalization matrix is given by $\boldsymbol{D}_n = \mathsf{diag} \left( n^{1/2}, n \right)$, $\boldsymbol{X}(s) = \left[ 1, B_2(s) \right]^{\prime}$.

Moreover, $d_2 (s)$ is a step-function such that, $d_2 (s) = 1$ if $s > 0.5$ and $0$ otherwise. Then, it holds that

and the variance under the null hypothesis is estimated to be

Testing for a Structural Break in Cointegrating Regression

Consider the following predictive regression model

where $y_t$ is a scalar vector and $\boldsymbol{x}_t \in \mathbb{R}^{p \times 1}$ is $p-$dimensional regressor vector. Furthermore, consider the case in which $\boldsymbol{\Phi}_T= \boldsymbol{I}_p$. Define the vector $\boldsymbol{u}_t = \big( \epsilon_t, \boldsymbol{\eta}_t^{\prime} \big)^{\prime}$ and the following long-run covariance matrices

such that

Similarly, we define

Furthermore, we define with $\boldsymbol{\Omega}_{ \epsilon \eta } = \boldsymbol{\Omega}_{ \epsilon \epsilon } - \boldsymbol{\Omega}_{ \epsilon \eta } \boldsymbol{\Omega}_{ \eta \eta }^{-1} \boldsymbol{\Omega}_{ \eta \epsilon }$ and $\boldsymbol{\Lambda}_{ \eta \epsilon }^{+} = \boldsymbol{\Lambda}_{ \eta \epsilon } - \boldsymbol{\Omega}_{ \eta \eta } \boldsymbol{\Lambda}_{ \eta \eta }^{-1} \boldsymbol{\Lambda}_{ \eta \epsilon }$. Then, we estimate the cointegrating regression model via OLS and obtain the residuals $\hat{\epsilon} = y_t - \boldsymbol{\beta}^{\prime} \boldsymbol{z}_t$, where $\boldsymbol{z}_t = \big( 1, \boldsymbol{x}_t^{\prime} \big)^{\prime}$ and $\boldsymbol{\beta} = \big( \beta_0, \boldsymbol{\beta}_1^{\prime} \big)$. Moreover, we define the following vector of residuals

Notice that $\Delta \bar{z}$ represents the sample mean of $\Delta z$. Furthermore, we denote with $\hat{u}_t$ the estimated vector of residuals as well as with $\hat{\boldsymbol{\Omega}}$ and $\hat{\boldsymbol{\Lambda}}$ the corresponding estimators of the long-run covariance matrices. We also define the transformed dependent variable as

Therefore, we obtain the fully modified (FM) estimator as below

The associated residual vector is given by $ \hat{ \epsilon }^{+} = y_t^{+} - \boldsymbol{\beta}^{+} \boldsymbol{x}_t$. Furthermore, we define the following

Define the following test statistic

The F-test is a test for a structural break at a known point which is asymptotically distributed under the null hypothesis as a $\chi^2-$distributed random variable with $p$ degree of freedoms. However, when the break-point is unknown, then the test is $F_{ \mathsf{sup} } = \mathsf{sup}_{ k \in \mathcal{B} } F_k$, where $\mathcal{B} \in (0,1)$. Moreover, a framework for structural break detection in predictive regression models is proposed by katsouris2023predictability, katsouris2023structural.

Unit root and Cointegrating limit theory

FM transformation in Cointegrating Regression

A challenging issue in the cointegrating regression literature is that when $u_t$ is uncorrelated with $v_t$ and hence uncorrelated with $x_t$, it follows that (i) $\lambda_{uv} = \boldsymbol{0}$, $\Delta_{uv} = \boldsymbol{0}$, and (ii) $B_u(r)$ is independent of $B_v(r)$. Furthermore, because of the independence between the Brownian motions $B_u(r)$ and $B_v(r)$ in this case, one can condition on $B_v(r)$ to show that the limiting distribution of $T \left( \widehat{\beta} - \beta \right)$ is a zero mean Gaussian mixture. Therefore, one can also show that the $t$ and Wald statistics for testing the hypotheses about $\beta$ have the usual $\mathcal{N}(0,1)$ and chi-square limits assuming serial correlation in $u_t$ is handled using consistent robust standard errors. On the other hand, when regressors are endogenous, the limiting distribution of $T \left( \widehat{\beta} - \beta \right)$ is obviously more complicated because of the correlation between $B_u(r)$ and $B_v(r)$ and the presence of nuisance parameters in the vector $\Delta_{vu}$.

Thus, we can no longer condition on $B_v(r)$ to obtain an asymptotic normal result and $\Delta_{vu}$ introduces an asymptotic bias. In particular, inference is difficult in this situation because nuisance parameters cannot be removed by simple scaling methods. Therefore, the FM-OLS estimator of phillips1990statistical (see, also phillips1990asymptotic) is employed to asymptotically remove $\Delta_{vu}$ and to deal with the correlation between $B_u(r)$ and $B_v(r)$. In practise FM-OLS estimation requires the choice of bandwidth and kernel. While bandwidth and kernel play no role asymptotically when considering the consistency results for $\widehat{\Omega}$ and $\widehat{\Delta}$, in finite samples they affect the sampling distributions of $\widehat{\theta}^{+}$ and thus of $t$ and Wald statistics based on the FM-OLS estimator of $\widehat{\theta}^{+}$. In other words, to obtain an approximation for $\widehat{\theta}^{+}$ that reflects the choice of bandwidth and kernel, the natural asymptotic theory to use is the fixed-b theory. However, fixed-b theory has primarily been developed for models with stationary regressors, which means that some additional work is required to obtain analogous results for cointegrating regressions. A major difference is that the first component of $\widehat{\eta}_t$, that is, $\widehat{u}_t$, is the residual from a cointegrating regression, which leads to dependence of the corresponding limit partial sum process on the number of integrated regressors and the deterministic components.

Specifically, the FM-OLS estimator for the cointegrating predictive regression model is defined as below

where $y_t^{+} = y_t - \hat{\boldsymbol{\Omega}}_{12} \hat{\boldsymbol{\Omega}}_{22}^{-1} D \boldsymbol{x}_t$. Then, the asymptotic distribution of the FM-OLS estimator

which is the same as the asymptotic distribution of the CCR estimator. Now, we may use the residuals

The FM-OLS is optimal in the sense that is has the same limiting distribution as Gaussian maximum likelihood. Then estimators can estimated with long-run covariance matrices.

The OLS residuals can be written as below

Consider the partial sum processes $S^{u}_t = \sum_{i=1}^t u_i$ and $S^{x}_t = \sum_{i=1}^t x_i$, which implies that

where

In particular, these results show that eliminating nuisance parameters is not easy unless $\boldsymbol{x}_t$ is strictly exogenous. Thus, the endogeneity aspect of the model implies that there is a dependence between the nonstationary regressor $X_t$ and the stationary error $u_t$.

FM estimation with nearly integrated regressors

In particular, the autoregressive root of the regressor which is parametrized with the local-to-unity specification allows to capture the near unit root behaviour of many predictor variables and is less restrictive than the pure unit root specification. Moreover, the OLS estimator of $\beta$ does not have an asymptotically mixed normal distribution due to the correlation between $u_t$ and $v_t$. From the paper of hjalmarsson2007fully the fully modified OLS estimator is given by

Define with $\omega_{11.2} = \omega_{11} - \omega_{21}^2 \omega_{22}^{-1}$ and $B_{1.2} = B_1 - \omega_{21} \omega_{22}^{-1} B_2 \equiv BM ( \omega_{11.2} )$. Then, as $T \to \infty$, the following limiting distribution holds

Furthermore, under the assumption that both $u_t$ and $v_t$ are martingale difference sequences it can be shown that the OLS estimation of the augmented regression

yields an estimator of $\beta$ with an asymptotic distribution identical to that of $\hat{\beta}^{+}$.

Dynamic Seemingly Unrelated Cointegrating Regression

We follow the framework proposed by mark2005dynamic. More specifically, we consider a fixed number of $N$ cointegrating regressions each with $T$ observations.

Moreover, we denote with $\boldsymbol{u}_{t} = \left( u_{1t}, ..., u_{NT} \right)^{\prime}$ and $\boldsymbol{e}_{t} = \left( e^{\prime}_{1t} ,... , e^{\prime}_{Nt} \right)$ and with $\boldsymbol{w}_t = \left( \boldsymbol{u}_{t}^{\prime}, \boldsymbol{e}_{t}^{\prime} \right)^{\prime}$ an $N(k+1)$dimensional vector with the orthonormal Wold moving average representation such that $\boldsymbol{w}_t = \boldsymbol{\Psi}(\boldsymbol{L}) \boldsymbol{\epsilon}_t$. Notice that $\boldsymbol{\epsilon}_t$ is a martingale difference sequence, such that $\mathbb{E} \left[ \boldsymbol{\epsilon}_t \right] = \boldsymbol{0}$ and $\mathbb{E} \big[ \boldsymbol{\epsilon}_t \boldsymbol{\epsilon}_t^{\prime} \big] = \boldsymbol{I}_k$ and finite fourth moments. Therefore, the endogeneity problem shows up as correlation between the $i-$th equilibrium error $u_{it}$ and potentially an infinite number of leads and lags of the first differenced regressors from all of the equations of the system such as $\Delta \boldsymbol{x}_{jt} = \boldsymbol{e}_{jt}$ for $(i, j) \in \left\{ 1,..., N \right\}$.

Consider the following formulation of the model

Let $\boldsymbol{y}_t = \left( y_{1t},..., y_{Nt} \right)^{\prime}$, $\boldsymbol{u}_{t} = \left( u_{1t},..., u_{Nt} \right)^{\prime}$ and $\boldsymbol{\beta} = \left( \boldsymbol{\beta}_1^{\prime} ,..., \boldsymbol{\beta}_N^{\prime} \right)^{\prime}$, $\boldsymbol{\delta} = \left( \boldsymbol{\delta}_1^{\prime} ,..., \boldsymbol{\delta}_N^{\prime} \right)^{\prime}$. Moreover we have the following matrices

Then, the equations can be stacked together in a system as below

Moreover, denote the long-run covariance matrix of $\boldsymbol{u}_t$ by $\boldsymbol{\Omega}_{uu}$. Then, an expression for the DSUR estimator with known $\boldsymbol{\Omega}_{uu}$ is given by

This, it can be shown that $\hat{ \boldsymbol{\beta} }_{DSUR}$ is asymptotically mixed normal.

Hypothesis testing for the linear restrictions $\boldsymbol{R} \boldsymbol{\beta} = \boldsymbol{r}$, where $\boldsymbol{R}$ is a $q \times Nk$ matrix of constants and $\boldsymbol{r}$ is a $q-$dimensional vector of constants. In practise, the testing hypothesis can be also written as below

Notice that the above formulation of the null hypothesis provides a way for conveniently formulating a test of homogeneity restrictions on the cointegrating vectors.

Then, the DSUR Wald test statistic is expressed as below

where

The Wald statistic $\mathcal{W}_{DSUR}$ is asymptotically distributed as a chi-square variate with $q$ degrees of freedom under the null hypothesis. Moreover, in applications we replace $\boldsymbol{\Omega}_{uu}$ with a consistent estimator, $\hat{ \boldsymbol{\Omega} }_{uu} \overset{ p }{ \to } \boldsymbol{\Omega}_{uu}$. Estimation of the long-run covariance matrix is discussed below. Such an estimator might be called a "feasible" DSUR estimator.

\paragraph{Two Step DSUR}

The first step purges endogeneity by least squares and the second step estimates $\boldsymbol{\beta}$ by running SUR on the least squares residuals obtained from the first-step regression. Under standard regularity conditions, the two-step DSUR estimator is asymptotically equivalent to the DSUR estimator $\hat{ \boldsymbol{\beta} }_{DSUR}$ discussed above. Stacking the equations together as $\hat{ \boldsymbol{y} }_t = \hat{ \boldsymbol{X} }_t^{\prime} \boldsymbol{\beta} + \hat{ \boldsymbol{u} }_t$ and running SUR gives the two-step DSUR estimator gives,

\paragraph{Restricted DSUR}

Next, we consider the estimation of the cointegration vector under the homogeneity restrictions $\boldsymbol{\beta}_1 = ... = \boldsymbol{\beta}_N = \boldsymbol{\beta}$. As in two-step DSUR, endogeneity can first be purged by regressing $y_{it}$ and each element of $\boldsymbol{x}_{it}$ on $\boldsymbol{z}_{pt}$. Let $\hat{y}_{it}$ and $\hat{x}_{it}$ denote the resulting regression errors. Therefore, the estimation problem becomes one of estimating $\boldsymbol{\beta}$, in the system of equations $\hat{y}_{it} = \hat{\boldsymbol{x} }_{it}^{ \prime } \boldsymbol{\beta} + \hat{u}_{it}$. Stacking the system equations together:

\paragraph{Asymptotic Properties} \

Let $\boldsymbol{W}(r)$ be a vector standard Brownian motion for $0 \leq r \leq 1$.

Notice that a functional central limit theory applies such that

with the long-run covariance matrix $\boldsymbol{\Omega} = \mathsf{diag} \big( \boldsymbol{\Omega}_{uu} , \boldsymbol{\Omega}_{ee} \big)$.

Due to the block diagonality of $\boldsymbol{\Omega}$, it can be easily verified that $\boldsymbol{B}_u$ and $\boldsymbol{B}_e$ are independent as discussed in the frameworks of mark2005dynamic and saikkonen1993estimation. Similarly phillips1999linear develop the asymptotic theory for regression models with nonstationary panel data for which the time series component is an integrated process and where both $T$ and $n$ are large. Recently, chen2023seemingly consider a framework for robust estimation of a Seemingly Unrelated Regression System for VAR models with explosive roots based on the approach of magdalinos2009limit and magdalinos2020least.

Predictive Regression Models

Regression with Integrated Regressors

Various examples related to integrated processes and the limit theory in regression with integrated predictors can be found in the book of banerjee1993co. Among others, related literature includes the papers of cavanagh1995inference, stock1991confidence, jansson2006optimal. Below we present an example (see, corresponding chapter from banerjee1993co), to shed light on the related asympotics.

The OLS estimator of $\beta$ is given by

Thus,

Note that we have a regression with an integrated regressor, since

Therefore, it follows that

In order to derive the limiting distribution of the model parameter in the case of integrated regressors, we shall first derive the limiting distribution of $\left( T^{-1} \sum_{t=1}^T x_t u_t \right)$.

To do this, we condition $u_t$ on $\epsilon_t$ as given below

Define $W_{\epsilon} (r)$ and $W_{v} (r)$ to be two independent Wiener processes on $\mathcal{C}[0,1]$. Therefore,

Substituting $x_t = x_{t-1} + \epsilon_t$ into the above gives

We have the following asymptotic results

Therefore, we have that

Furthermore, phillips1988asymptotic proved the following result

Therefore, under the null hypothesis $H_0: \beta = 0$,

Thus, the $t-$statistic, denoted as $\mathcal{T}_{\beta = 0}$ for testing the null hypothesis, $H_0: \beta = 0$,

has the following limiting distribution

The above limiting distribution indicates that the $t-$ratio of $\hat{\beta}$ does not follow a standard normal distribution unless $\phi = 0$, which in that cases implies that $x_t$ is exogenous for the estimation of $\beta$. In particular, when $\phi \neq 0$ then the first term of the above limiting distribution gives rise to second-order or endogeneity bias, which although asymptotically negligible in estimating $\beta$ due to super consistency, can be important in finite samples.

Optimal Inference in Predictive Regressions with Persistent Regressors

Following jansson2006optimal consider the case where $\left\{ ( y_t, x_t^{\prime} \right\}$ is generated by the predictive regression model with the local to unity specification for the autocorrelation matrix.

Let $\hat{\Omega}$ be a consistent estimator of

which is the long-run variance of $\left( \epsilon_{yt}, \Psi(L) \epsilon_{xt} \right)^{\prime}$.

Therefore, within our setting we aim to compare the efficiency of the likelihood ratio test and the Wald type statistics as structural break detectors for the predictive regression model with persistent regressors at an unknown break-point. In order to study the asymptotic behaviour of the individual components we consider limit theorems involving invariance principles for partial sums within the moderate deviations framework developed by phillips2007limit.

Notice that the stochastic process $x_{ \floor{ns}, n} := x_{ \floor{ns} } / d_n$ on the Skorohod space $\mathcal{D}[0,1]$ converges weakly to a Gaussian process $G(s)$. Furthermore, on a suitably expanded probability space there exists a process $\left( x_{ t, n}^{0}, 1 \leq t \leq n \right) =_d \left( x_{ t, n}, 1 \leq t \leq n \right)$ such that

Moderate Deviations from a Unit Root in Autoregressive Processes

It is possible to approximate the partial sum process on $\mathcal{D} [ 0, + \infty )$ for $x_{[ n^{\gamma} \bullet ]}$ by the Stieltjes integral

Therefore for each $\gamma \in (0,1)$ and $c < 0$, we have that

\color{black}

Notice that with the seminal paper of phillips2007limit, the authors introduce limit theorems and invariance principles for moderate deviations from unit root. The particular framework is employed in time series models such as autoregressive and predictive regression models with predictors assumed to be generated with the LUR specification. A strong approximation over the interval $[ 0, n^{1 - \alpha} ]$ for the partial sum process of i.i.d errors can be constructed via an expanded probability space with a Brownian motion $W(.)$ with variance $\sigma^2$ for which

Therefore, we have that one of the main results of the paper which provide a uniform approximation of $n^{- \alpha / 2} y_{ \floor{ n^{\alpha} . } }$ by $V_{ n^{\alpha } }$ on $[ 0, n^{1 - \alpha} ]$. For each $\alpha \in (0,1)$ and $c < 0$ we have that

The importance of the above result lies in the fact that an embedding of the random element $V_{ n^{\alpha } } (t)$ to the linear diffusion $J_c(t) := \displaystyle \int_0^t e^{ c(t-r) } dB(r)$ is possible. Using integration by parts it can be shown that,

Moreover, by Lemma 3.1 we obtain the following expression

Therefore, we obtain that

We note that the limit theory is established through a combination of a functional law to a diffusion and a central limit law to a Gaussian random variable. A more immediate consequence is the limit law of the random element $y_{ \floor{ n^{\alpha} . } }$ on the original probability space. For all $j = 0,..., \floor{ n^{1 - \alpha} } - 1$ and $p \in [0,1]$ we obtain $\frac{1}{ n^{ \alpha / 2} } y_{ \floor{ n^{\alpha } j} + \floor{ n^{\alpha } p} } $.

We analyze each of the two terms separately. The term containing the block components can be written:

Taking the inner integral along $[0, r] = [ 0, \floor{ n^{1 - \alpha} } ]$ we have, up to $o_p(1)$ that

where the remainder term $R_n$ is shown in the Appendix to be $o_p(1)$.

The second integral on the right side can be evaluated directly to obtain

Furthermore, we obtain that

Then, $\Psi_n^{\prime} \to_{\text{a.s}} \Psi = \sum_{ j=1 }^{ \infty } \rho^{-j} u_j^{ \prime }$, and it follows by the Shorokhod representation theorem that $\Psi_n \to_d \Psi$. Then, by joint weak convergence of $\Psi_n$ and $Z_n$ it follows that $\left( Z_n, \Psi_n \right) \Rightarrow \left( Z, \Psi \right)$ as $n \to \infty$, with $Z =_d \Psi$. The limiting random variables $\Psi$ and $Z$ can be shown to be independent by modifying Anderson's argument adjusted for weakly dependent errors. The idea is that, as $n \to \infty$, $Z_n$ can be approximated by the first $\floor{ L_n }$ elements of the sum $\sum_{ j = 1}^n \rho^{ -j } u_j$ whereas $\Psi_n$ can be approximated by the last $\floor{ L_n }$ elements of the sum $\sum_{ j = 1}^n \rho^{ - (n-j) -1} u_j$, where $\left( L_n \right)_{ n \in \mathbb{N} }$ is a sequence increasing to $\infty$ with $L_n \leq n / 3$ for each $n$.

Accordingly, we define

We may further approximate

We have that, $\Psi_n - \Psi_n^{*} = \sum_{ k = \floor{ L_n } }^n \rho^{ -k } u_{ n - k + 1}$, and so

Similar to the sample variance, the asymptotic behaviour of the sample covariance is partly determined by elements of the time series $y_{t-1} u_t$ that do not belong to the block components

Obtaining limits for the block components and the remaining time series separately in a method similar to that used for the sample variance will work. For the limit as $\alpha \to 1$, we have that $n^{1 - \alpha} \to 1$, and so $\floor{ n^{1 - \alpha }} = 1$ for $\alpha = 1$, in which case $j = 0$ in the aforementioned blocking scheme. Therefore, the invariance principle proposed by phillips1987time $n^{ - 1 / 2} y_{ \floor{np} } \Rightarrow J_c(p)$ on $\mathcal{D} \left( [0,1] \right)$ which yield the usual local to unity limit result

In other words continuity in the limit theory cannot be achieved at the boundary with the conventional local to unity asympotics, at least without using the blocking construction.

\paragraph{Supplementary Results: Main Proofs}

An immediate consequence of the above Proposition is the following identity. For each $x \in [ 0, n^{1 - \alpha} ]$ and $m \in \mathbb{N}$

Since, by Proposition A3 and the fact that $\sum_{t=1}^n \left| \rho_n \right|^t = \mathcal{O} \left( n^{\alpha} \right)$,

\paragraph{Proof of Lemma 4.1}

Consider the stochastic process

is a stationary process with autocovariance function given by

Then by Theorem of Inrangimov and Linnik (1971)

Therefore,

provided that $\sum_{ h = 0}^{ \infty } \left| \gamma_{\xi} (h) \right| < \infty$.

To show summability of the covariance function $\gamma_{ \xi }$, write

Once the fitted stationary process is obtained, the variance $\Sigma$ of $\left\{ w_t \right\}$ can be consistently estimated as

\paragraph{Summary} Although it is not shown explicitly, the idea of stationary transformations can be used also for the essentially singular model to construct the regression which yields asympotically invariant, with respect to nuisance parameters, tests.

The use of moderate deviation principles are useful in understanding the various asymptotic theory results that apply across the parameter space of the autoregressive coefficient of the nonstationary autoregressive model. In particular, for the case of quantile autoregressions and quantile predictive regression models a relevant framework is proposed by katsouris2022asymptotic.

Multivariate Predictive Regression Model

Consider the multivariate predictive regression model with lag-augmentation

where $\left( \boldsymbol{y}_{t}, \boldsymbol{x}_t \right) \in \mathbb{R}^{p \times 1 }$ are $p-$dimensional vectors for $t \in \left\{ 1,..., n \right\}$ and $\left\{ \boldsymbol{A} \right\} \in \mathbb{R}^{p \times p }$ are $(p \times p)$ matrices of coefficients and $\boldsymbol{C} = \mathsf{diag} \left\{ c_1,..., c_p \right\}$ is a $(p \times p)$ diagonal matrix with the coefficients of persistence. The degree of persistence of the regressors is determined by the value of the unknown persistence coefficients $c_i$'s which are assumed to be positive constants. Moreover, the common exponent rate for the persistence coefficients take values, $\alpha = 1$, which covers the local-unit-root (LUR) regressors and $\alpha \in (0,1)$ the case of mildly integrated (MI) regressors as defined by magdalinos2009limit and kostakis2015Robust. When $\alpha = 0$ we assume the presence of stationary regressors, for a suitably restricted coefficient matrix $\left( \boldsymbol{I}_p - \boldsymbol{C}_p \right)$. A more recent application is given by yang2020testing.

Moreover, it holds that $B(r) = \omega W(r)$ for some $r \in [0,1]$ and $\omega := \sum_{j = - \infty}^{\infty} \gamma_u(j) = C(1)^2 \sigma^2$, the one-sided covariance of $v_t$. Moreover, by definition it holds that $W_c(r) = \int_0^r e^{ c(r-s)} dW(s)$. Notice that, expression (ref) provides an invariance principle for the partial sum process of $\boldsymbol{w}_t$, which implies that the partial sum of $\boldsymbol{w}_t$, weakly converges to the stochastic quantity $\boldsymbol{\Sigma}_{w}^{1/2} \boldsymbol{W}(r)$, that is, a $p-$dimensional Brownian process with covariance matrix $\boldsymbol{\Sigma}_{w}$. Within our context we employ the local-unit-root limit theory as proposed by the seminal paper of phillips1987time, such that $\boldsymbol{x}_{[nr]} / \sqrt{n} \Rightarrow \boldsymbol{J}_c(r)$, where $\boldsymbol{J}_c(r)$ is a $p-$dimensional Gaussian process defined as below

that satisfies the Black-Scholes differential equation $d \boldsymbol{J}_c(r) \equiv \boldsymbol{C}_p \boldsymbol{K}_c(r) + d \boldsymbol{B}_v(r)$, with $\boldsymbol{K}_c(r) =0$, implying also that $\boldsymbol{K}_c(r) \equiv \sigma_v \boldsymbol{J}_c(r)$, where $\displaystyle \boldsymbol{J}_c(r) = \int_0^r e^{(r-s)\boldsymbol{C}_p} d \boldsymbol{W}(s)$. Moreover, since in this paper we consider time series nonstationarity in the form of regressors exhibiting mildly integradeness, then we also employ the asymptotic results proposed by phillips2007limit since the invariance principle $x_{[nr]}$ requires appropriate normalization. Notice that $\boldsymbol{K}_c(r)$ represents the Ornstein-Uhlenbeck, (OU) process, which encompasses the unit root case such that $\boldsymbol{J}_c(r) \equiv \boldsymbol{B}_v(r)$, for $\boldsymbol{C}_p = 0$. Moreover, the assumption of a local-unit-root specification for the predictors of the model such that $\boldsymbol{x}_t = \left( \boldsymbol{I}_p - \frac{ \boldsymbol{C}_p }{n^{ \alpha } } \right)$ $\boldsymbol{x}_{t-1} + \boldsymbol{v}_t$ and more specifically by allowing the autocorrelation coefficient to be of the form $\boldsymbol{\Phi}_n = \left( \boldsymbol{I}_p - \frac{ \boldsymbol{C}_p }{ n^{ \alpha } } \right)$, permits to consider the $\boldsymbol{K}_c(r)$ Gaussian process given by (ref), as a stochastic approximation.

Proof. Expanding the expression $x_t z_t^{\prime}$, vectorizing and summing over $t \in (1,...,n)$

We have that

Therefore, for $\delta < \gamma$, which holds since $\gamma = 1$ and $\delta \in (0,1)$, then we obtain

Note for example, that

Thus,

Part (L2). Since, $J_c(r) = \int_0^r e^{-c(r-s)} dB_v(s)$ and using expression (14) in PM for the univariate predictive regression with single regressor we have that

The corresponding result for the the univariate predictive regression with multiple regressors is as below

Recall that by definition the functional $\underline{J}_c(r) = \int_0^r e^{(r-s)C} d \underline{B}(s)$, is considered to be a vector diffusion process and satisfies the stochastic differential equation system

Therefore, using the initial condition as well we can also write

Since $\underline{J}_c(r)$ is a Gaussian process and for fixed $r$ then it can be easily proved that the finite dimensional distribution (see e.g., phillips1988regression)

where $\Omega$ the $p \times p$ covariance matrix of $\underline{B}(r)$ the $p-$dimensional Brownian motion.

Part (L3). The weakly convergence result to a mixed Gaussian distribution below

shows that the limit distribution of $T^{ - (1 + \delta )/2} \sum_{t=1}^T \left( z_{t-1} \otimes u_t \right)$ is Gaussian with mean zero and covariance matrix equal to the probability limit of $T^{ - (1 + \delta )/2} \sum_{t=1}^T \left( z_{t-1} \otimes u_t \right)$, which is equal to $V_{zz} \otimes \Sigma_{ uu }$, where $V_{zz} := \int_0^{\infty} e^{r C_z} \Omega_{vv} e^{r C_z} dr$. Similarly the limit distribution of $T^{ - (1 + \delta )/2} \sum_{t=1}^T \left( x_{t-1} \otimes u_t \right) \Rightarrow N \left( 0 , V_{xx} \otimes \Sigma_{ uu } \right)$, where $V_{xx} := \int_0^{\infty} e^{r C} \Omega_{vv} e^{r C} dr$, is proved in Lemma 3.3 of magdalinos2009limit.

We prove the weakly convergence result given above for the univariate predictive regression (applied to both the case of single versus multiple regressors). See also case 2 of Lemma B3 in Appendix of KMS, where the result for the multivariate predictive regression is proved. We consider the partial sum process $S_t = \sum_{j = 1}^t z_j$ where $z_j$ is the IVX instrument. Therefore, we obtain $\Delta S_{t-1} \equiv S_{t-1} - S_{t-2} = z_{t-1}$, i.e.,

The first term of the above expression is considered to be asymptotically negligible, i.e., $o_p(1)$. Thus, we consider the asymptotic behaviour of the second term by expanding the expression further. We have that $\displaystyle \Delta x_t = \frac{C}{T} x_{t-1} + v_t$, implying

Therefore, the above gives the following

We prove the limiting behaviour of the above terms using the convergence of the martingale expression (16) in PM and the fact that $x_t$ is a local to unity process.

Then, the following weakly convergence of (ref) and (ref) follows

Predictive Regression with Multiple Regressors

Consider the predictive regression model such as

where $y_t$ is a scalar vector and $\boldsymbol{x}_{t-1}$ is a $k-$dimensional vector of regressors. Denote with $\boldsymbol{\eta}_t = \big[ u_t, \boldsymbol{v}_t^{\prime} \big]^{\prime}$.

Then, the predictive regression above includes a model intercept and furthermore we do not demean the random variables in order to remove the model intercept. In addition, we assume that $\boldsymbol{R}_T \neq \boldsymbol{I}_T$ which would correspond to the integrated regressor case.

Fully Modified OLS Under the aforementioned conditions the functional law

holds for partial sums of $\boldsymbol{\eta}_t$. Define the partition $\boldsymbol{B} = \big( B_u, \boldsymbol{B}_v^{\prime} \big)^{\prime}$.

When $\Omega_{uu.v}$ has full rank, the rate of convergence of the FM-OLS estimator is determined by the rates of weak convergence of the sample covariances and the rate of nonparametric estimation of $\Omega$ and $\Gamma^{+}$ does not apply any role. In particular, the OLS estimator $\widehat{A}^{+} = Y^{\prime} X \left( X^{\prime} X \right)^{-1}$ and employs corrections for endogeneity in the regressor $x_t$, leading to the transformed dependent variable

and a bias correction term involving

which is constructed in the usual way using consistent nonparametric estimators of submatrices of the long run and one sided long run quantities $\Omega$ and $\Gamma^{+}$. In particular, the FM-OLS estimators removes asymptotic bias and increases efficiency by correcting both the long run serial correlation in $u_t$ and endogeneity in $x_t$ causes by the long run correlation between $u_{0t}$ and $u_{xt}$.

Testing for Structural Break in Predictive Regression

In this section, we present the framework of georgiev2018testing. The approach of the particular paper is to consider between two important test statistics which are commonly used in detecting structural change, that is, the sup Wald test of andrews1993tests and the Cramer-von-Mises type statistics of nyblom1989testing. On the one hand the sup Wald test proposed by andrews1993tests is motivated by alternatives where the parameters display a small number of breaks at deterministic points in the sample, while under the alternative hypothesis of the Cramer-von-Mises test statistic the coefficients of the model are random and slowly evolve through time. Using a fixed regressor wild bootstrap procedure their framework allows for both conditional and unconditional heteroscedasticity in the data.

Consider the time-varying predictive regression model

where $\epsilon_{yt}$ is a mean zero innovation process and $x_t$ is an observed process, generated by

where $\rho_x := \left( 1 - \frac{c_x}{T} \right)$, with $c_x \leq 0$ which allows to consider two types of persistence properties for the regressors of the model, that is, (i) strongly persistence unit root, and (ii) local-to-unit root.

\paragraph{Stochastic Coefficient Variation}

The assumption of time-varying parameters in the predictive regression model implies that the time-varying components of the model parameters, that is, $s_{\alpha t}$ and $s_{\beta t}$ follow (near-) unit root processes expressed as below

where the autocorrelation coefficient is expressed as below

which are unit-root or local-to-unit root autoregressive processes.

The non-stochastic coefficient variation mechanism includes the testing framework proposed by andrews1993tests. Under the assumption of one-time breakpoint, then $s_{ \alpha t}$ and $s_{ \beta t}$ are modelled as below

where we define $D_t \left( \floor{ \pi_0 T } \right) := \mathbf{1} \left( t \leq \floor{ \pi T } \right)$ with $\floor{ \pi T }$ denoting a generic shift point with associated break fraction $\pi$. We take the true shift fraction $\pi_0 \in \Pi$, where $\Pi = [ \pi_1, \pi_2]$ with $0 < \pi_1 < \pi_2 < 1$.

Consider the testing hypothesis given by

Moreover, consider the fitted regression model

Let $\hat{\sigma}^2 := T^{-1} \sum_{t=1}^T \hat{e}_t^2$, with $\hat{e}_t^2$ the OLS residual from the fitted regression given by (ref). Moreover, we denote with $\underline{X} = [ 1 \ x_{t-1} ]^{\prime}$.

Then, the LM statistic is constructed as

Note that the corresponding single parameter LM statistics are given by

where LM$_1$ corresponds to the test statistic relating the the intercept alone and LM$_2$ to the slope coefficient alone.

\paragraph{Supplementary Examples}

In this section we provide additional examples related to testing for parameter instability in predictive regression models. Further derivations can be found in the studies of katsouris2023predictability, katsouris2023structural, katsouris2023testing. First, consider the case in which the predictive regression model includes a model intercept but is assumed to be stable before conducting any statistical inference on the existence of a structural break on the slopes of the model.

Main Aspects on Predictive Regressions

A research question of interest is whether future values of one series $\left\{ Y_t \right\}$ can be predicted from lagged values of another series $\left\{ X_t \right\}$, where the two series obey the following model:

It is well known that the OLS estimate of $\beta$ is biased when the errors $\left\{ u_t, v_t \right\}$ are correlated, with the amount of bias increasing as $\rho$ get closer to the unit boundary. There is currently no known theoretical justification that inference (e.g., a t-test) based on the bias-corrected OLS estimate will have improved size properties relative to the test based on the uncorrected OLS estimate. In practise, there are several studies in the literature that provide examples where even the use of the exact bias correction does not result in accurate finite-sample inference because the resulting t-statistic can still be very far from normal. Therefore, from an inferential point of view, issues such as the bias of point estimates may be irrelevant in small samples, and thus it is also crucial to examine the likelihood. At first glance, it may seem surprising that a likelihood ratio test (LRT) may provide well-behaved hypothesis tests in situations when the t-statistic does not because the to test statistics are closely related (see, phillips2014restricted).

In this section, we consider the estimation of the predictive regression model from the statistical point of view, motivating the use of the restricted likelihood, which is free of the nuisance intercept parameter and hence able to imitate the likelihood of the no-intercept univariate model with its attendant small curvature. A curvature-related approach to tackle the predictive regression problem for the model was also taken by jansson2006optimal. Due to the presence of model intercepts $( \mu, \eta )$, we are motivated to seek a likelihood that does not involve the location parameters and yet possesses small curvature properties similar to those of the model with known location parameters. Thus, the restricted likelihood has exactly these properties. In practise, the idea of restricted likelihood was originally proposed by KS (1970) precisely as a means of eliminating the effect of nuisance location parameters when estimating the parameters of the model covariance structure in a linear model.

We focus our attention on the use of the restricted likelihood for carrying inference on $\beta$ in the bivariate predictive regression model. Notice that in the particular framework, the restricted likelihood is the exact likelihood of the vector $\big\{ ( Y_t - Y_{t-1} )_{t=2}^n, ( X_t - X_{t-1} )_{t=2}^n \big\}$. In practise, the restricted maximum likelihood (REML) estimates are found to be asymptotically efficient. A standard method for testing the two-sided hypothesis, $H_0: \beta = 0$ versus $H_A: \beta \neq 0$ is the LRT, which compares the log-likelihood evaluated at the unrestricted estimates of the parameters to the log-likelihood evaluated at the parameter estimates obtained under the restriction that the null hypothesis $H_0: \beta = 0$ is true. Then, it is easy to generalize the REML likelihood in two directions that are both of practical interest. One generalization is to the case where the predictor series is a multivariate AR(1) process. The other generalization is to the case where the univariate predictor follows a higher order AR process.

In other words, the inclusion of an intercept in the model causes the likelihood ratio to lose this property, thus pointing to the intercept as the source of the problem. Therefore, this motivates the use of the restricted likelihood, which is free of the nuisance intercept parameter and hence able to imitate the likelihood of the no-intercept univariate model with its attendant small curvature. Thus, we are indeed able to obtain theoretical results that demonstrate the the LRT based on the restricted likelihood (RLRT) has good finite-sample performance for both estimation and inference in this context.

\paragraph{Multivariate Regression}

Assume that the data $\big( Y_1,..., Y_n, \boldsymbol{X}_0,..., \boldsymbol{X}_n^{\prime} \big)$ follow

where $u_t = \boldsymbol{\phi}^{\prime} \boldsymbol{v}_t + e_t$, such that $\big( e_t, \boldsymbol{v}_t^{\prime} \big)^{\prime} \sim \mathcal{N} \big( \boldsymbol{0}, \mathsf{diag} \left( \sigma_e^2, \boldsymbol{\Sigma}_v \right) \big)$ is an i.i.d series and $\boldsymbol{R}$ is a $(k \times k)$ matrix with all eigenvalues less than unity in absolute value. Define with $\boldsymbol{\Sigma}_v \equiv \mathsf{Var} ( \boldsymbol{v}_t )$ and $\boldsymbol{\Sigma}_{ \boldsymbol{X} } \equiv \mathsf{Var} \left( \boldsymbol{X}_{t} \right)$ as

and define with

Then, the REML log-likelihood up to an additive constant for the model is given by

where

and

Furthermore, since the dimension of the parameter space is very large in the vector case, it is not feasible to obtain a result such as Theorem 3 in the most general case. However, in the case where $\boldsymbol{R}$ is a diagonal matrix and where $\big( \sigma_e^2, \boldsymbol{\phi}, \boldsymbol{\Sigma}_v \big)$ are assumed known with $\boldsymbol{\Sigma}_v$ diagonal, we are able to obtain the following result on the finite-sample behaviour of the RLRT for testing $\mathbb{H}_0: \boldsymbol{\beta} = \boldsymbol{0}$.

\paragraph{Higher Order Autoregressive Regressors}

Let the observed data $\big( Y_1,..., Y_n, X_{-p+2},..., X_n \big)$ follow

where $u_t = \phi v_t + e_t$ and $\big( e_t , v_t \big) \sim \mathcal{N} \big( \boldsymbol{0}, \mathsf{diag} \left( \sigma_e^2, \sigma_v^2 \right) \big)$ are an $\textit{i.i.d}$ series. Moreover, assume that all the roots of the polynomial $z^p - \sum_{s=1}^p z^{p-s} \alpha_s$ lie within the unit circle. Define with

where $X_i = \frac{1}{n} \boldsymbol{1} \boldsymbol{X}_i$.

Resampling for Time Series Regressions

Seminal studies for aspects related to implementation of resampling methods in time series analysis include politis1994stationary, paparoditis2001tapered, paparoditis2003residual among others.

Preliminary Theory on Bootstrap Methodologies

Let $\left\{ X_n \right\}_{ n \in \mathbb{Z}}$ be a stationary sequence of random variables with common continuous distribution function $F(t) = \mathbb{P} \left( X_0 \leq t \right)$ on a probability space. Assume $0 \leq X_0 \leq 1$. The empirical process is defined as

where $F_n(t)$ is the empirical distribution of $\left\{ X_n \right\}$. We have that $B_n(t) \to_d B(t)$ in the spaces $D[0,1]$ endowed with Skorohod topology where $B$ is a Gaussian process specified by $\mathbb{E} \left( B(t) \right) = 0$ for every $t$, and for every $t$ and $s$,

In addition,

The block-based bootstrap estimators of the mean and empirical process are defined as below. Let $k$ and $l$ be two integers such that $n = kl$. Let $T_{n1}, T_{n2},..., T_{nk}$ be i.i.d random variables each having uniform distribution on $\left\{ 1,...,n \right\}$. Define the triangular array $\left\{ X_{ni}, 1 \leq i \leq n + l \right\}$by $X_{ni} = X_i$ for $1 \leq i \leq n$ and $X_{ni} = X_{i - n }$ for $1 \leq i \leq n + l$. In other words, we extend our sample of size $n$ by another $l$ observations, namely, $X_1,...,X_l$. Then, the bootstrapped estimator of the mean is defined as below

and the bootstrapped estimator of the empirical process is given by

Then, the boostrapped empirical process is defined as

We denote with $\mathbb{E}^{*}$, $var^{*}$ the moments under the conditional probability measure $\mathcal{P}^{*}$ induced by the resampling mechanism, that is, $\mathcal{P}^{*}$, is the conditional probability given $\left( X_1,...,X_n \right)$. In this paper, we focus on the weak convergence of $B_n^{*}(t)$ to $B(t)$ in the Skorohod topology on $D[0,1]$, almost surely.

\paragraph{Simulating the bootstrap distribution}

We follow the steps below to obtain the approximation to the bootstrap distribution $G_n \left( x, \widehat{F}_n \right)$. We consider that $b = 1,...,B$ are the bootstrap replications.

Now, how the sample in Step 1 of Algorithm 1 is drawn, depends on the form of the estimator $\widehat{F}_n$. Two main bootstrap methodologies are commonly used:

Bootstrap in regression models

Consider the following regression model for $i = 1,...,n$

\paragraph{Pair bootstrap} The pairs bootstrap is a direct extension of the iid bootstrap and is based on the bivariate EDF $\widehat{F}_n (x,y)$ for the sample $(X_1, Y_1),..., (X_n, Y_n)$ are a random sample with bivariate cdf F. The pairs bootstrap builds the bootstrap sample $( X_1^{*}, Y_1^{*} ),...,( X_n^{*}, Y_n^{*} )$ by drawing pairs with replacement from the sample $(X_1, Y_1),..., (X_n, Y_n)$. It is crucial that $(X_i, Y_i)$ are kept together as a pair, otherwise if $\textbf{Y}^{*}$ and $\textbf{X}^{*}$ were drawn separately the bootstrap would assume there is no relation between the two.

On Bootstrap Asymptotics

The bootstrap is a method for estimating the distribution of an estimator or test statistic by resampling one's data or a model estimated from the data. The methods that are available for implementing the bootstrap and the improvements in accuracy that it achieves relative to first-order asymptotic approximations depend on whether the data are an i.i.d random sample or a time series. If the data are iid, the bootstrap can be implemented by sampling the data randomly with replacement or by sampling a parametric model of the distribution of the data. The distribution of a statistic is estimated by its empirical distribution under sampling from the data or parametric model (see, hardle2003bootstrap). Therefore, the situation is more complicated when the data are a time series because bootstrap sampling must be carried out in a way that suitably captures the dependence structure of the DGP.

Why the bootstrap provides Asymptotic Refinements

The term asymptotic refinements refers to approximations to distribution functions, coverage probabilities that are more accurate than those of first-order asymptotic distribution theory (see, beran1988prepivoting). We focus on the distribution function of the asymptotically $N(0,1)$ statistic $T_n = n^{1 / 2} \left( \theta_n - \theta \right) / s_n$.

Let $\widehat{ \mathcal{P} }$ denote the probability measure induced by bootstrap sampling, and let $\widehat{T}_n$ denote a bootstrap analogue of $T_n$. If the data are i.i.d, then it suffices to let $\widehat{ \mathcal{P} }$ be the empirical distribution of the data. Bootstrap samples are drawn by sampling the data $\left\{ X_i : i = 1,...,n \right\}$ randomly with replacement. If $\left\{ \widehat{X}_i: i = 1,...,n \right\}$ is such a sample, then

where $\hat{\theta}_n = \theta( \hat{m}_n )$ and $\hat{m}_n = n^{-1} \sum_{i=1}^n \widehat{X}_i$, and $\widehat{s}_n$ is obtained by replacing the $\left\{ X_i \right\}$ with $\left\{ \hat{X}_i \right\}$ in the formula for $s_n^2$.

The Block Bootstrap