The modified conditional sum-of-squares estimator for fractionally integrated models

The modified conditional sum-of-squares estimator for fractionally integrated models

\fancyhf \fancyhead[L]{ \sf This is the Manuscript Accepted by the Journal of Econometrics. \newline \textcopyright \ 2026. This manuscript version is made available under the \href{https://creativecommons.org/licenses/by-nc-nd/4.0/}{CC-BY-NC-ND 4.0 license}.} \thispagestyle{fancy}

\pagestyle{plain}

Introduction

Fractionally integrated autoregressive moving average (ARFIMA) models are applied in a wide range of fields for describing long-memory phenomena, witness inter alia the economic and political as well as the natural sciences; see hassler2019time and hualde2021frac for general treatments. One particular variant of this model class that has recently gained popularity is the so-called type-II ARFIMA model, which truncates the fractional integration operator and allows both stationary and non-stationary processes to be described, see for example nielsen2004efficient, robinson2005distance and johansen2008representation. A popular choice for estimating this model is the conditional sum-of-squares (CSS) estimator whose main appealing features are that it is computationally straightforward and that the memory parameter can be estimated consistently as long as it lies in an arbitrary compact interval on the real line. It was introduced by li1986fractional in the context of stationary fractionally integrated models. Subsequent papers allowed for non-stationary models, see for instance beran1995maximum and velasco2000whittle. Local consistency proofs were provided by tanaka1999nonstationary, nielsen2004efficient and robinson2006conditional. Global consistency was proved by hualde2011gaussian and nielsen2015asymptotics in a model without deterministic components. Only recently, hualde2020truncated,hualde2021truncated derived global consistency and the asymptotic normality of the CSS estimator in a model with deterministic components, such as a constant or a trending term. Empirical applications include hualde2011gaussian for aggregate income and consumption data and johansen2016role for opinion poll data.

While the literature dealing with asymptotic inferences in the context of parametric ARFIMA models is well-developed, some issues still require attention. One such concern pertains to the small sample performance of the CSS estimator. Despite the widespread use of the CSS estimator little is currently known about the impact deterministic terms have on the properties of the estimator of the memory parameter in small samples. Early on, chung1993small and cheung1994maximum conducted simulation studies and found that the inclusion of a constant term in the model can substantially increase the small-sample bias and mean squared error (MSE) of the estimated memory parameter. lieberman2005expansions and johansen2016role are among the few theoretical contributions to shed light on the issue. lieberman2005expansions derive the Edgeworth expansion of the memory parameter for the Gaussian maximum likelihood estimator in stationary fractional time series model. johansen2016role investigate the impact of observed and unobserved initial values on the bias of the memory parameter estimator in a non-stationary fractional time series model. Neither paper, however, includes short-run dynamics in its model. In addition, we are not aware of any related work that simultaneously tackles both stationary and non-stationary processes.

The purpose of the present paper is therefore to add to this literature and analyse the small-sample bias of the CSS estimator in a type-II fractional model with short-run dynamics and constant term from an analytical, empirical and simulation point of view. In particular, our analysis reveals that incorporating the level parameter into the model introduces an additional bias in the CSS estimator. This bias is due to a biased score which is particularly pronounced when the data is stationary. We will suggest what we call the modified conditional sum-of-squares (MCSS) estimator which (i) is easy to compute, (ii) removes the leading bias term and (iii) allows much more accurate small-sample inference.

To do so, we will interpret the constant term as nuisance parameter and draw on a large literature on bias correction. laskar1998modified provide an overview of this literature. We build on the approach to dealing with nuisance parameters initiated by conniffe1987expected and mccullagh1990simple and recently applied by bartolucci2016modified and by martellosio2020adjusted, i.e.\ we adjust the score function so that its expectation equals zero. The idea is as follows: We find a stochastic higher-order expansion of the estimator as a function of the derivatives of the profile likelihood, cf.\ johansen2016role,lawley1956general. The expansion is simplified by approximating the derivatives by their leading terms. This allows the expectation of the estimator to be evaluated explicitly. We show that a leading component of the bias is attributable to the nonzero expectation of the score. By premultiplying the objective function by a suitable modification term then results in the expected score evaluated at the true parameter to be equal to zero, thereby mitigating the bias of the estimator.

It is important to emphasise that this paper tackles the correction of the CSS estimator's bias that is due to the estimation of the unknown constant term in the model, henceforth referred to as “unknown-level bias”. Of course, there may be other, additional, sources of bias, yet they are not the focus of the present investigation. In particular, as we will discuss below, the estimator is also subject to what we call “intrinsic bias”, i.e.\ the bias also present if the constant term in the model is known. Moreover, what we will refer to as “misspecification bias” is due to conditioning on unobserved pre-sample values, as considered by johansen2016role and hualde2020truncated. While there is no way around the intrinsic bias in our model, we abstract initially from the misspecification bias for expositional clarity by making a suitable simplifying assumption. Subsequently, our results are extended to a more general model setup that includes unobserved pre-sample values.

The main contributions of this paper to the literature are threefold: First, we examine our MCSS estimator in type-II ARFIMA($p_1$,$d$,$p_2$) models with constant term and compare it to the standard CSS estimator. In particular, we derive its exact bias and we show that it is consistent and asymptotically normally distributed. The results generate new insights into the sources of bias and into bias correction of other models nested in our setup, such as purely fractional ARFIMA models and stationary ARMA($p_1$,$p_2$) models. Secondly, we re-visit three classical datasets that have in the past been described by ARFIMA($p_1$,$d$,$p_2$) models with constant term, namely the post-second World War real GNP data, the extended Nelson-Plosser dataset, and the Nile data, by applying our MCSS estimator to estimate the long-memory parameter and the short-run dynamics. All three time series are short and therefore warrant the use of small-sample bias corrections. Our conclusion sheds new light on the interpretation of these datasets. Thirdly, this paper paves the way to extending the analysis of small-sample bias from univariate type-II ARFIMA modes to panel settings, see also the contributions of robinson2015efficient and schumann2023role.

The remainder of the paper is organised as follows. In Section (ref) we present the MCSS estimator for a parametric fractional time series model. In Section (ref) we conduct a simulation study to examine the small sample properties of the estimators. Section (ref) presents the empirical illustrations. Section (ref) contains concluding remarks. All proofs are relegated to the appendix.

The modified conditional sum-of-squares estimator

The model

Consider a so-called type-II fractional process $z_t$, $t = 0,\pm 1,\pm 2,\ldots$, generated by the model

with $\epsilon_t \sim \textit{IID}(0,\sigma^2)$, where $0 < \sigma^2 < \infty $. Here, $\Delta = 1 - L$ and $L$ are the difference and lag operators, respectively, and $d$ can take any value in $\mathbb{R}$. For any series $v_t$, real number $\zeta$ and time index $t \geq 1$, the so-called truncation operator $\Delta_+^{\zeta}$ is defined by $\Delta_+^{\zeta} v_t = \Delta^{\zeta} \{ v_t I(t \geq 1) \} = \sum_{i = 0}^{t-1} \pi_{i}(-\zeta) v_{t-i},$ with $I(\cdot)$ being the indicator function and with $\pi_{i}(a)$ denoting the coefficients in the usual binomial expansion $\Delta^{-a} = \sum_{i = 0}^{\infty} \pi_i(a) L^i$, where $\pi_{0}(a) = 1$ and $\pi_{i}(a) = (i - 1 + a) \pi_{i-1}(a) / i,$ for $i = 1, 2, \ldots$, see e.g.\ hassler2019time. The parameter $d$ in (ref) is known as the memory parameter or the fractional parameter. A consequence of the indicator function in the definition of the truncation operator is that $z_t = 0$ for all $t \leq 0$. The process $z_{t}$ has been widely applied in the literature, see marinucci2000weak,marinucci2001semiparametric, robinson2003cointegration, nielsen2004efficient, shimotsu2005exact, robinson2005distance and johansen2008representation, among others. An extension of this setup is introduced by johansen2016role who allow for $N_0$ unobserved pre-sample values of $z_t$ by defining $ \Delta_{-N_0}^{\zeta} v_t = \sum_{i = 0}^{t+N_0-1} \pi_i(-\zeta) v_{t-i}$. As a result, (ref) would become $z_t = \Delta_{-N_0}^{-d} \epsilon_t$. Setting $N_0 = 0$, the original truncation operator is recovered; the greater the value of $N_0$, the longer the burn-in period before the process is observed.

Two comments on the memory parameter are of interest: First, its range is commonly divided into a “stationary” and a “non-stationary” region: $d < 1/2$ and $d \geq 1/2$, respectively. Since the definition of the truncation operator implies that $z_t = 0$ for $t \leq 0$ the process $z_t$ is in fact not covariance stationary when $d < 1/2$ and $d \neq 0$. However, it may be considered asymptotically stationary for any such $d$. To see this, consider the so-called type-I fractional process $\tilde{z}_t = \Delta^{-d} \epsilon_t$ which is known to be covariance stationary for any $d < 1/2$. marinucci1999alternative observe that for $|d| < 1/2$, $ E\left( z_t - \tilde{z}_t \right)^2 = O(t^{2d -1}) $ as $t \rightarrow \infty$, and hence the difference between $\tilde{z}_t $ to $z_t$ vanishes. Although marinucci1999alternative consider only $|d| < 1/2$, their result actually holds for any $d < 1/2$. This follows from Stirling's approximation and johansen2016role. This asymptotic equivalence prompts us to retain the terminological dichotomy between stationarity and non-stationarity. Secondly, it is worth noting that even for $d \geq 1/2$, i.e.\ in the non-stationary region, the truncation operator ensures that the process $z_{t}$ is well-defined in the mean-square sense, see johansen2008representation and hualde2011gaussian.

While the model in (ref) covers a wide range of dynamics, it seems unsuitable for many empirical applications for two key reasons. First, it implies that $E(z_t) = 0$. Secondly, the simple IID errors assumed in model (ref) are overly restrictive. In order to make our model more widely applicable, we therefore complement the model in (ref) by a constant term $\mu$ and replace $\epsilon_t$ by $u_t$ to yield

where $\omega$ is a lag polynomial and captures the short-run dependence structure parametrically, given by

with $\varphi$ being an unknown $p\times1$ vector, $\omega_0(\varphi) = 1$, $|\omega(s;\varphi)| \neq 0$ for $|s| \leq 1$, and $\sum_{j = 0}^{\infty} |\omega_j(\varphi)| < \infty$. More precise conditions on $\omega$ will be specified below. The representation of $u_t$ in (ref) as a MA$(\infty)$ model is common in the literature and considered by, among others, hualde2011gaussian and hualde2020truncated,hualde2021truncated. As a consequence of the constant term, $E(x_t) = \mu$ for $t \geq 1$. Note that equation (ref) implies that $x_t = 0$ for $t \leq 0$ although equation (ref) allows for the process $u_t$ to have started in the infinite past. Indeed, the literature has so far considered several different specifications of pre-sample values, as we will discuss further below.

One well-known special case of $u_t$ in (ref) is an ARMA$(p_1,p_2)$ model which is given by

where $\beta(L;\varphi)$ is the AR polynomial of order $p_1$ and $\alpha(L;\varphi)$ is the MA polynomial of order $p_2$. It is assumed that the polynomials do not have common roots and that their roots lie outside the unit circle. Then (ref), (ref) and (ref) is an ARFIMA$(p_1,d,p_2)$ model. Another special case of $u_t$ in (ref) is bloomfield1973exponential's (1973) exponential spectrum model, see robinson1994efficient and hassler2019time.

Following hualde2020truncated, we make the following assumptions on the model's short-run dynamics and on the admissible parameter space. We use the notation $\vartheta = (d,\varphi')'$. True parameter values are denoted by the subscript 0, i.e.\ $\vartheta_0 = (d_0,\varphi_0')'$ and $\mu_0$.

Assumption (ref) states that the errors $\epsilon_t$ are IID, but it does not restrict them to normality. In fact, the IID assumption can be weakened to martingale difference series as in hualde2020truncated,hualde2021truncated but for the sake of convenience we keep this condition simple. Assumption (ref) covers the admissable parameter spaces and the true parameter values. Note that we exclude the boundary case of $d_0 = 1/2$. The reason is that our analysis of the CSS estimator's bias takes its consistency for granted. To the best of our knowledge, hualde2020truncated,hualde2021truncated are the only papers to show consistency in a model with constant term, yet they exclude the value of $d_0 = 1/2$.\footnote{Specifically, hualde2020truncated,hualde2021truncated study CSS estimation in models with generalised polynomial or power-law trends whose coefficients and exponents are unknown, and they exclude the boundary case $d_0 = 1/2$ for technical reasons. A constant term is, of course, a known deterministic structure, yet proving consistency of the CSS estimator in this model is beyond the scope of our analysis.} Assumption (ref)$(i)$-$(iv)$, which ensures the identification of the short-term dynamics, is standard in the literature on parametric short-memory models since its introduction by hannan1973asymptotic. Assumption (ref)$(v)$-$(vii)$ serve as additional regulatory conditions necessary to establish the asymptotic distribution theory. We refer to the papers by hualde2011gaussian, nielsen2015asymptotics, hualde2020truncated,hualde2021truncated for a detailed discussion of Assumption (ref). Importantly, it is satisfied for the stationary and invertible ARMA model and also the exponential spectrum model of bloomfield1973exponential.

For the purpose of expositional clarity we also make the following assumption in our initial analysis. Subsequently, our main results will be generalised to a model setting in which this restriction is ignored.

Assumption (ref) boils down to not only the obervations $x_t$ in (ref) but also the unobserved error terms $u_t$ and $\varepsilon_t$ in (ref) being zero before $t = 1$. Assumption (ref) reflects the definition of the CSS estimator in that the latter conditions on all pre-sample information. Importantly, it is not required for the consistency or asymptotic normality of our MCSS estimator, nor does it have an impact on the correction of the unknown-level bias. Instead, as will become plain in the generalisation presented in Section (ref), it helps to disentangle the unknown-level bias from the model misspecification bias, i.e.\ from the bias that is due to the CSS estimator not making use of pre-sample information. Several papers in the literature have employed an assumption similar to Assumption (ref), see e.g.\ SibbertsenEtAl18 and robinson2020estimation. Alternative initialisation schemes also exist: the paper by johansen2016role allows for a finite number of unobserved pre-sample values, while hualde2011gaussian and hualde2020truncated assume that both $u_t$ and $\varepsilon_t$ have started in the infinite past. Indeed, this latter generalisation will be used in Section (ref) below to dispense with Assumption (ref) and to extend the main results of our analysis.

The conditional sum-of-squares estimator

We now discuss the conditional sum-of-squares (CSS) estimator of the parameters in model (ref)-(ref). This is the estimator considered by e.g.\ hualde2011gaussian who, however, look at a model without the constant term. We distinguish the case in which $\mu$ is unknown from that in which it is known. As will be seen in Section (ref) below, the CSS estimator may also be motivated as a conditional maximum likelihood estimator under the assumption of Gaussian innovation terms $\epsilon_t$, as in johansen2016role and hualde2020truncated.

Consider a sample of observations for $t = 1, \ldots, T$. For any $(\vartheta,\mu)$ in the admissible parameter space, define the residuals $\epsilon_t(\vartheta,\mu) = \phi(L;\varphi) \Delta_+^{d}(x_t-\mu)$. The CSS objective function is then given by

where we define the convoluted coefficient

with $ \kappa_{0t}(d) = \Delta_{+}^{d} I(t \geq 1)= \sum_{n = 0}^{t-1} \pi_n(-d) = \pi_{t-1}(1-d), $ cf.\ johansen2016role.

Since $L(\vartheta,\mu)$ in (ref) is quadratic in $\mu$ we can concentrate it. Unsurprisingly, the CSS estimator of $\mu$ for fixed $\vartheta$ is given by

Substituting $\hat \mu (\vartheta)$ into (ref) yields the profile (or concentrated) CSS function

Note that we use asterisks to emphasise that we are dealing with a profile objective function. The resulting CSS estimator of $\vartheta = (d,\varphi')'$ is given by

As discussed in Section (ref), the model effectively conditions on $x_t = 0$, for $t \leq 0$. For this reason, hualde2011gaussian and hualde2020truncated prefer to call the estimator in (ref) the truncated sum-of-squares estimator. hualde2020truncated show that if $x_t$ is generated by (ref)-(ref) and if Assumptions (ref) to (ref) hold, then, as $T \rightarrow \infty$, $\hat{\vartheta} \overset{P}{\rightarrow} \vartheta_0$ and $ \sqrt{T} (\hat{\vartheta} - \vartheta_0 ) \xrightarrow{D} N(0_{p+1},A^{-1}), $ where $A$ is given in (ref). Note that the authors exclude the boundary case of $d_0 = 1/2$ from their analysis, see also the discussion of our Assumption (ref) above. We take this consistency of the CSS estimator as basis for the following analysis of its bias.

A few remarks about the estimator $\hat{\mu}(\vartheta)$ in (ref) are instructive. For $\vartheta = \vartheta_0$ we have that $\hat{\mu}(\vartheta_0)-\mu_0 = \sum_{t = 1}^T \epsilon_t c_t(\vartheta_0) / \sum_{t = 1}^T c^2_t(\vartheta_0) ,$ which has mean zero and variance $\sigma_0^2 (\sum_{t = 1}^T c^2_t(\vartheta_0))^{-1}$. In the stationary region, i.e.\ when $d_0 < 1/2$, this variance goes to zero because then $\sum_{t = 1}^T c^2_t(\vartheta_0)$ diverges in $T$, see Lemma (ref). As opposed to that, in the non-stationary region, i.e.\ when $d_0 > 1/2$, this variance does not converge to zero because then $\sum_{t = 1}^T c^2_t(\vartheta_0)$ is bounded in $T$, see Lemma (ref). This is the reason why $\hat{\mu}(\hat{\vartheta}) \xrightarrow{P} \mu_0$ only if $d_0 < 1/2$, see hualde2020truncated for the proof.

For comparison, we also analyse the situation where the true $\mu_0$ is known. The CSS estimator for this model can be derived by substituting $\mu_0$ into (ref) to have

such that

This estimator is considered by hualde2011gaussian and nielsen2015asymptotics who show that if $x_t$ is generated by (ref)-(ref) and if Assumptions (ref) to (ref) hold, then $\hat \vartheta_{\mu_0}$ is consistent, too, attaining the same limiting distribution as $\hat{\vartheta}$. In other words, the distribution does not depend on whether $\mu$ is known or needs to be estimated.

The modified profile likelihood

A central concern in this paper is to investigate the bias of $\hat{\vartheta}$ in (ref) and of $\hat{\vartheta}_{\mu_0}$ in (ref). This will be done in Section (ref) below. It will turn out that the expectation of the CSS estimators is a function of the expectation of the score functions, or first derivatives, of $L^*(\vartheta)$ and $L_{\mu_0}^*(\vartheta)$ evaluated at $\vartheta = \vartheta_0$, respectively. The present section therefore examines the bias of the two scores and builds on an approach by mccullagh1990simple to correct for it.

To that end, it will be instructive to interpret the CSS objective in (ref) as a log-likelihood function, as do johansen2016role and hualde2020truncated. Assuming for the moment that $\epsilon_t \sim \textit{NID}(0,\sigma^2)$, the Gaussian log-likelihood of $x_t$ in (ref), conditional on $x_t$ = $0$ for $t\leq 0$, is given by $l (\vartheta,\mu,\sigma^2) = - T \log (\sigma^2 ) /2 - \sum_{t = 1}^T ( \phi(L;\varphi) \Delta_+^{d} (x_t-\mu ) )^2 / (2\sigma^2). $ Throughout the paper, we omit additive constants in log-likelihood functions for notational simplicity. Maximising $l (\vartheta,\mu,\sigma^2)$ with respect to $\sigma^2$ yields $ \hat{\sigma}^2(\vartheta,\mu) = T^{-1} \sum_{t = 1}^T ( \phi(L;\varphi) \Delta_+^{d} (x_t-\mu ) )^2 $ and the resulting profile log-likelihood is

Maximising $\ell (\vartheta,\mu)$ further with respect to $\mu$ results in $\hat \mu (\vartheta)$ in (ref) and the profile log-likelihood function becomes

Clearly, the estimator of $\vartheta$ resulting from maximising (ref) is identical to that obtained by minimising (ref) since

So, the CSS objective $L^*(\vartheta)$ can be seen as a negative non-logged profile likelihood even though we do not impose Normality on the error term $\epsilon_t$ in Assumption (ref). As the maximum likelihood estimator of $\vartheta$ is asymptotically efficient, see hualde2020truncated, so is the CSS estimator $\hat \vartheta$ in (ref). The same can of course be said of $\hat \vartheta_{\mu_0}$ in (ref) since the profile CSS objective $ L^*_{\mu_0}(d)$ in (ref) can be obtained from (ref) by replacing $\mu$ by its known value $\mu_0$ such that $ L^*_{\mu_0}(\vartheta) = T \exp ( - 2 \ell(\vartheta,\mu_0)/T ) /2 . $

We will in the present section therefore interpret $L^*(\vartheta)$ in (ref) as a profile likelihood. As such, it is not a genuine likelihood, for it is not directly based on observable quantities, see barndorff1983formula and severini2000likelihood. Instead, it is a function of the maximum likelihood estimators of $\mu$ and $\sigma^2$ which are treated as if they were the true parameter values. In large samples, the concentration procedure has relatively minor effects, yet chung1993small showed in Monte Carlo simulations that in small samples it leads to a strong bias in $\hat{\vartheta}$. This is because profile likelihoods do not necessarily possess the same properties as genuine likelihoods. It is well-known that, under classical regularity conditions and with a fixed number of regressors, the score of the profile likelihood is biased. In particular, its expectation is $O(1)$, see kalbfleisch1973marginal, mccullagh1990simple and liang1995inference. The following theorem derives the bias of the scores of $L^*(\vartheta)$ and $L^*_{\mu_0}(\vartheta)$. The proof will be given in Appendix (ref) and a generalisation that dispenses with Assumption (ref) in Theorem (ref) of Appendix (ref). We use the notation $D_i f(\vartheta) = D_i f(d, \varphi)$ to denote the first derivative of a function $f (\vartheta) = f(d, \varphi)$ with respect to parameter $i \in \{\vartheta, d, \varphi\}$.

Clearly, the scores $\mathit{D_dL}^*(d_0,\varphi_0)$ and $\mathit{D_{\varphi}L}^*(d_0,\varphi_0)$ are biased. In addition, the expectation of the score in (ref), i.e.\ the score with respect to $d$, is not uniform in $d_0$. For $d_0 < 1/2$ it diverges at the rate of $\log(T)$, while it is $O(1)$ for $d_0 > 1/2$. This is not true for the expectation of the score in (ref), i.e.\ the score with respect to $\varphi$, which is $O(1)$ uniformly in $d_0$. The rationale behind this is that the score bias measures the relative strengths of the level parameter and the stochastic component. The score bias with respect to $d$ gauges the strength of the level parameter relative to the fractional dynamics, whereas the score bias with respect to $\varphi$ evaluates the strength of the level parameter in relation to short-run dynamics. In the non-stationary region, i.e.\ when $d_0 > 1/2$, we recall that $\mu$ is not consistently estimated, see the discussion in Section (ref). The reason is that the stochastic component $\Delta_+^{-d} u_t$ in (ref) dominates the deterministic component $\mu$. Hence, the bias in the score is less influenced by $\hat{\mu}(\hat{\vartheta})$, resulting in the expected scores being $O(1)$ for such $d_0$. On the other hand, if $d_0 < 1/2$, $\mu$ is consistently estimated and $\hat{\mu}(\hat{\vartheta})$ plays a more important role in the bias of the scores and especially for the score with respect to the fractional dynamics. This is reflected by the expected score in (ref) being $O(\log(T))$ for such $d_0$. Recall that, at the beginning of this section, we mentioned that the bias of the CSS estimator is a function of the score bias. One might be tempted to think, from (ref) and (ref), that in the stationary region the bias of $\hat{d}$ will be of a larger order of magnitude than the bias of $\hat{\varphi}$. Yet this turns out not to be true. As will be shown below, the biases of $\hat{d}$ and $\hat{\varphi}$ are functions of not only their own score biases but, instead, of a weighted sum of both score biases. This will lead to the order of the bias of the short-run dynamics to be the same as that of the memory parameter. The situation for $L_{\mu_0}^*(d,\varphi)$ in (ref) is somewhat different. Although, technically speaking, $L_{\mu_0}^*(d,\varphi)$ is also a profile likelihood due to the substitution of $\hat{\sigma}^2$ for $\sigma^2$, its scores are unbiased, as shown in (ref) and (ref).

This discussion highlights the need for a modification of the profile likelihood function such that it behaves more like a genuine likelihood in terms of score unbiasedness. This modification will eliminate the bias of the CSS estimator $\hat \vartheta$ stemming from the presence of the unknown nuisance parameter $\mu$, as will be seen in Section (ref). The idea of modifying the profile likelihood to obtain score unbiasedness is in fact not new and was previously discussed by mccullagh1990simple. martellosio2020adjusted, for instance, implement this idea for a spatial model.

To develop the idea of an unbiased score we follow mccullagh1990simple by considering for simplicity a profile log-likelihood in a scalar parameter $\vartheta$. Recentering its score yields, say, $ D \ell^*_{a}(\vartheta) = D \ell^*(\vartheta) - a(\vartheta), $ where $ \ell^*(\vartheta)$ denotes, as before, the profile log-likelihood function and where $a(\vartheta)$ is an adjustment function only depending on $\vartheta$. Then mccullagh1990simple require that $ E \left( D \ell^*_{a}(\vartheta_0) \right) = 0, $ which implies that

for all $\vartheta_0$. Finally, they call $\ell^*_{a}(\vartheta) = \int_{\Theta} D \ell^*_{a}(t) dt$ the adjusted profile log-likelihood, which is subsequently maximised w.r.t.\ $\vartheta$.

Note that mccullagh1990simple further adjust $D\ell^*_{a}(\vartheta)$ to make it information unbiased, i.e.\ by making its variance equal the negative expectation of the derivative of the score. While these adjustments may improve the efficiency of the estimator, they are not addressed in the present paper because they do not affect the location of the zeros of $D\ell^*_{a}(\vartheta)$.

The adjustment function $a(\vartheta)$ can in principle be computed from (ref). Yet this calculation is challenging as can be seen by rewriting (ref) as

In special cases, the computation of $a(\vartheta)$ may indeed be feasible, see for instance the spatial model with Gaussian errors considered by martellosio2020adjusted. In general, however, the evaluation of (ref) is not straightforward. In the sequel, we therefore present an argument of how the expectation of the fraction can be circumvented. To that end, we consider the profile CSS objective function $L^*(\vartheta)$ as basis for the adjustment. Indeed, mccullagh1990simple in their Remark 3 allude to the possibility of using an objective function other than the profile log-likelihood $\ell^*(\vartheta)$ for deriving an adjustment.

With a view to framing the approach of mccullagh1990simple in terms of $L^*(\vartheta)$, note that the adjusted profile log-likelihood can be written as

with $A(\vartheta) = \int_{\Theta} a(t) dt$. Based on the relationship between $ \ell^*(\vartheta)$ and $L^*(\vartheta)$ in (ref), we can write (ref) as

Clearly, maximising the adjusted profile log-likelihood $ \ell^*_a (\vartheta)$ in (ref) is equivalent to minimising the adjusted profile CSS objective $L^*_{a}(\vartheta) = T \exp ( -2 \ell^*_{a}(\vartheta)/ T) / 2 $ which, upon substituting in $ \ell^*_a (\vartheta)$, is

It is important to note that while the adjustment in (ref) is additive, it is multiplicative in (ref). Recall that $a(\vartheta)$ in (ref), and thus $A(\vartheta)$ in (ref), is difficult to compute. We therefore define, as an alternative, the modified profile CSS objective function

where the multiplicative modification term $m(\vartheta) > 0$ depends only on $\vartheta$. The corresponding score function is the first derivative of (ref):

As do mccullagh1990simple, we now require that our objective function is score unbiased, i.e.\ that the score function in (ref) satisfies $E \left( \mathit{DL}^*_{m}(\vartheta_0) \right) = 0. $ Using (ref) and the fact that $ D \log (m(\vartheta)) = D m(\vartheta) /m(\vartheta)$ as well as $m(\vartheta) > 0$, it follows that this condition is equivalent to

It will be seen in the next section that the evaluation of (ref) is straightforward, as opposed to the evaluation of (ref).

The modified conditional sum-of-squares estimator

Re-write now the condition in (ref) in terms of a parameter vector $\vartheta$ so $D_{\vartheta}L^*(\vartheta)$ denotes the vector of derivatives of $L^* (\vartheta)$ w.r.t.\ $\vartheta$. The condition can now be used for finding the modification term $m(\vartheta)$ for the modified profile CSS objective in (ref): First, it is shown in Lemma (ref) that the expectation of $D_{\vartheta}L^*(\vartheta_0)$ equals

where $c_{t}(\vartheta_0)$ is given in (ref). It is also shown in Lemma (ref) that $E\left( L^*(\vartheta_0) \right) = \sigma^2_0 (T-1)/2$. Consequently, from (ref), we have

Integrating and exponentiating (ref) yields $ m(\vartheta) = e^k ( \sum_{t = 1}^T c^2_t(\vartheta) )^{1/\left(T-1 \right)}, $ where $k$ is the constant of integration.

The modified profile CSS objective function in (ref) is thus given by the product of $ m(\vartheta)$ and $ L^*(\vartheta)$. Clearly the argument that minimises the likelihood is invariant to the choice of $k$ so that, in the sequel, we can set $k = 0$ without loss of generality. The resulting modified profile CSS objective function is then given by

We call the argument that minimises $L_m^*(\vartheta)$ the modified conditional sum-of-squares (MCSS) estimator and denote it by $\hat{\vartheta}_{m}$, i.e.\,

The modification term, with $k = 0$ imposed, is

Two important properties of the modification term $m(\vartheta)$ are stated in the following lemma. See Appendix (ref) for the proof.

The property in (ref) implies that the modification term $m(\vartheta)$ acts as a penalty in the minimisation of the modified profile likelihood $L_m^*(\vartheta)$ through inflating $L^*(\vartheta)$ by $m(\vartheta)$. The property in (ref) ensures that $m(\vartheta) \rightarrow 1$ such that the asymptotic properties of the MCSS estimator $\hat \vartheta_m$ are the same as those of the CSS estimator $\hat \vartheta$ in (ref). This is desirable because the CSS estimator is efficient under Gaussianity, as argued in Section (ref). The asymptotic properties of $\hat \vartheta_m$ are summarised for completeness in the following theorem and are proved in Appendix (ref). Note that no use is made of Assumption (ref).

The intuitive explanation of Theorem (ref) follows from noticing that

where use was made of the definition $L^*_{m}(\vartheta_0)$ in (ref) and the asymptotic behaviour of $m(\vartheta)$ in (ref) of Lemma (ref). Since $ L^*(\vartheta_0)$ in (ref) is $O_P(T)$, the second summands have no influence on the asymptotic distribution of $\hat{\vartheta}_{m}$. For the bias, however, the latter terms require further analysis, which is carried out below in Section (ref).

We now provide an illustration of how the modification term behaves for some selected parameter values. For simplicity, we consider a purely fractional model, i.e.\ we set $\omega(L;\varphi) = 1$ in model (ref)-(ref) so that, effectively, $\vartheta = d$. The corresponding modification term $m(d)$ is then given by (ref) with $\phi(L;\varphi) = 1$. This modification term is plotted in panel (a) of Figure (ref) for some illustrative values of $d$ and $T$. Four important observations can be made: First, recall from (ref) that the modification term $m(d)$ penalises the CSS objective $L^*(d)$ through inflating it by the factor $m(d)$. It appears from the plot that in the stationary region, i.e.\ when $d < 1/2$, $m(d)$ inflates $L^*(d)$ more than in the non-stationary region, i.e.\ when $d \geq 1/2$. This is a reflection of the fact that the order of the bias in the score is larger in the stationary region, as was argued in (ref) of Theorem (ref). Secondly, it is plain that when $d = 1$ the bias caused by estimating the constant term $\mu$ is the smallest, as predicted in Lemma (ref). Thirdly, even for a moderately large sample of size $T = 256$, $m(d)$ still turns out to be substantial in the stationary region, implying that the corresponding bias in the score is large. Fourthly, the negative slope of the modification term $m(d)$ for $d < 1$ implies that the minimum of $L^*_m(d)$ is shifted to the right of that of $L^*(d)$. This is illustrated in panel (b) of Figure (ref) which displays a Monte Carlo simulation of the CSS and MCSS objective functions. The DGP is stationary and corresponds to the model in (ref)-(ref) with $\epsilon_t \sim \textit{NID}(0,1)$, $\omega(L;\varphi_0) = 1$, $d_0 = 0.2$ and $\mu_0 = 0$. The sample size is $T = 64$ and the number of replications is 10,000. On display is the Monte Carlo average of the simulated $L^*(d)$, $L^*_{\mu_0}(d)$ and $L_{m}^*(d)$. The solid line represents the Monte Carlo average of $L^*(d)$: it can be seen that the CSS estimator underestimates the true $d_0 = 0.2$ on average. The dash-dotted line represents the Monte Carlo average of $L_{\mu_0}^*(d)$, which takes the constant term as known. This estimator is, on average, close to $d_0$. The dotted line represents the Monte Carlo average of $L_{m}^*(d)$, whose minimum is shifted to the right of that of $L^*(d)$. It therefore corrects for the distortion in $L^*(d)$ caused by estimating $\mu$.

Relationship with other modifications

There is a large literature on correcting the bias of maximum likelihood caused by the presence of unknown nuisance parameters. Seminal contributions include barndorff1983formula who proposed the modified likelihood function, and cox1987parameter who contributed the idea of the conditional profile likelihood by approximating the modified likelihood function. Both modifications result in modified profile likelihoods that are approximately score unbiased, see liang1987estimating and cox1994inference. It is therefore illuminating to investigate how our MCSS objective, with an expected score exactly equal to zero, relates to alternative approaches to bias-reduction, or how our modification term $m(\vartheta)$ compares to alternative adjustments. This section discusses two such ideas.

First, reconsider the adjusted profile log-likelihood $\ell^*_a (\vartheta)$ proposed by mccullagh1990simple and derived in Section (ref). Denote the corresponding estimator by $ \hat{\vartheta}_{a} = \operatorname*{argmax}_{\vartheta \in \Theta} \ell^*_{a}(\vartheta)$. The following corollary establishes the condition under which the adjusted profile log-likelihood estimator $\hat \vartheta_a$ is identical to our MCSS estimator $\hat \vartheta_m$. The proof follows from (ref) and (ref) and is omitted.

However, there is no particular reason to expect (ref) to hold in general. Therefore, our modified estimator $\hat \vartheta_m$ will in general be different from the adjusted estimator $\hat{\vartheta}_{a}$ suggested by mccullagh1990simple.

Secondly, a modification term closely related to $m(\vartheta)$ in (ref) is the one discussed in an1993cox who implement the idea of cox1987parameter to adjust the log-likelihood function. The setup in an1993cox is different from ours, however: They consider a stationary Gaussian type-I ARFIMA$(p_1,d,p_2)$ process $\tilde{x}_t$ generated by the model $\tilde{x}_t = \mu + \Delta^{-d} u_t$ where $u_t$ is defined in (ref) and (ref) with $\epsilon_t\sim \textit{NID}(0,\sigma^2)$ and $|d| < 1/2$. This contrasts to our type-II process whose $d$ is also allowed to lie in the non-stationary region and whose error term $\epsilon_t$ is not assumed to be Normally distributed.

It will prove helpful to phrase the approach by an1993cox in matrix notation: Define the $T \times 1$ vector $\tilde{x} = (\tilde{x}_1,\ldots,\tilde{x}_T)'$ such that $\tilde{x} \sim N(\mu \iota, \sigma^2 \Sigma(\vartheta))$ where $\iota$ is a $T\times1$ vector of ones and $ \sigma^2 \Sigma(\vartheta)$ is the $T \times T$ variance-covariance matrix of $\tilde{x}$, see for instance hoskingdiff1981 for the elements of $\Sigma(\vartheta)$. The log-likelihood function is then given by $\tilde{\ell}(\vartheta,\mu,\sigma^2) = - \log( |\Sigma(\vartheta)| ) /2 - T \log( \sigma^2)/2 - ( \tilde{x} - \iota \mu )' \Sigma(\vartheta)^{-1} ( \tilde{x} - \iota \mu ) / (2 \sigma^2). $ Substituting in the maximum likelihood estimators $\hat{\mu}(\vartheta)$ and $\hat{\sigma}^2(\vartheta,\mu)$ yields the profile log-likelihood function

The modified profile log-likelihood function that an1993cox find is

It can be shown that the orders of magnitude of the last three summands on the right-hand side of (ref) are $O_P(1)$, $o(1)$ and $ O(\log(T))$, respectively. For the proof of the second summand, we refer to dahlhaus1989efficient. As for the third summand, denote by $f_{\vartheta}(\omega)$ the spectral density of $\tilde{x}_t$ such that $T^{-1} \log |\Sigma(\vartheta)| \rightarrow \frac{1}{2 \pi} \int_{-\pi}^{\pi} \log( f_{\vartheta}(\omega)/\sigma^2) d\omega $ which is equal to 0 by virtue of the well-known Szeg\H o--Kolmogorov formula, see Chan2006. The order of magnitude of the fourth summand follows from $\iota' \Sigma(\vartheta)^{-1} \iota = O(T^{1-2d+\varepsilon})$ for each $\varepsilon > 0$, cf.\ adenstedt1974large. The leading of the three summands is therefore the last one, its order of magnitude being $O(\log(T))$.

In order to compare $\tilde{\ell}^*_m (\vartheta)$ to our MCSS function $L_m^*(\vartheta)$ in (ref) we transform the latter again in a fashion similar to that in (ref) and define $ \ell^*_{m}(\vartheta) = - T \log ( 2 L_{m}^*(\vartheta) /T ) /2, $ such that $\ell^*_{m}(\vartheta) = \ell^*(\vartheta)- T \log ( m(\vartheta) ) /2, $ with $\ell^* (\vartheta)$ given in (ref). Using the definition of $m(\vartheta)$ in (ref) yields

since $T/(T-1) = 1 + 1/(T-1)$. As is shown in Lemma (ref), the second summand is of order $O (\log T)$.

Two observations are now instructive. First, the leading modification term in $\ell^*_m (\vartheta)$ is of the same order of magnitude as that in $\tilde \ell^*_m$ in (ref), namely $O(\log(T))$. Second, we note that the Cholesky factor of $\Sigma(\vartheta)^{-1}$ in (ref) is the GLS transformation matrix that filters out the correlation structure of the type-I error term $\Delta^{-d} u_t$. Similarly, in our setting, $\phi(L;\varphi)\Delta_{+}^{d}$ filters out the correlation structure of the type-II error $\Delta_+^{-d} u_t$. Indeed, if it were possible to replace $\Sigma(\vartheta)^{-1/2} \iota$ in (ref) by $\phi(L;\varphi)\Delta_{+}^{d} \iota$ we would obtain

using the definition of $c_{t}(\vartheta)$ in (ref). Let us emphasise again, however, that the approach by an1993cox, although asymptotically equivalent to ours, is based on a model that assumes stationary and Normally distributed data. In addition, it necessitates the computation of the $T \times T$ variance-covariance matrix, or its Cholesky factor, which is often onerous computationally.

Asymptotic biases

This section investigates the asymptotic biases of the estimators of $\vartheta$. Two questions are of central interest. First, by how much does the MCSS estimator $\hat{\vartheta}_m$ reduce the bias of the CSS estimator $\hat{\vartheta}$? Second, is the bias of the MCSS estimator $\hat \vartheta_m$ comparable to that of the CSS estimator with known $\mu_0$? To address both questions, we proceed in a similar fashion as do johansen2016role, involving two steps: first, we find asymptotic expansions of the estimators and, secondly, we approximate these expansions. Since there would be a risk of overloading the exposition with bulky notation, the treatment is here mainly of verbal-descriptive type. Technical details are contained in Appendix (ref). Additional intuition can also be gained from the discussion of the asymptotic biases in the simple model without short-run dynamics, presented in Section (ref).

With a view to deriving the asymptotic bias of the CSS estimator $\hat{\vartheta}$, take a second-order Taylor series expansion of $D_\vartheta L^* ( \hat \vartheta) = 0$ around $\vartheta_0$, which results in an expression involving the first three derivatives of $L^* ( \vartheta)$, denoted by $D_\vartheta L^* ( \vartheta)$, $D_{\vartheta\vartheta'} L^* (\vartheta) $, and $D_{\vartheta_i \vartheta \vartheta'} L^* (\vartheta) $ for $i = 1,\ldots, p+1$. Since the representation of the third derivative in particular is unwieldy, we refer the reader to equation (ref) for an explicit formula. We demonstrate in Lemmata (ref) and (ref) that the derivatives satisfy $D_\vartheta L^* ( \vartheta_0) = O_P(T^{1/2})$, $D_{\vartheta'\vartheta} L^* (\vartheta_0) = O_P(T)$, and $D_{\vartheta_i \vartheta \vartheta'} L^* (\vartheta_0) = O_P(T)$ for $i = 1,\ldots, p+1$. Importantly, the orders of magnitude hold uniformly in $d_0$, allowing us to treat the stationary and non-stationary region jointly. Using the Taylor expansion of $D_\vartheta L^* ( \hat \vartheta) = 0$ and the expressions of the derivatives of $L^* ( \vartheta)$ allows us to find $G_{1T}$ and $G_{2T}$ in $\hat{\vartheta}-\vartheta_0 = T^{-1/2} G_{1T} + T^{-1} G_{2T} + O_p(T^{-3/2})$ and, thence,

where $S_T(\vartheta_0 ) = - A^{-1} T^{-1} [ \sigma^{-2}_0 E \left( \mathit{D_{\vartheta}L}^*(\vartheta_0) \right) ]$ and $B_T (\varphi_0)$ are given in equations (ref) and (ref), respectively. The purpose of the decomposition in (ref) is to separate the expectation of the first derivative of $L^* (\vartheta)$ from the correlations between the derivatives: the former enters $S_T (\vartheta_0)$ only, while the latter are gathered in $B_T (\varphi_0)$. We therefore call $S_T (\vartheta_0)$ the score bias of $\hat \vartheta$ and $B_T(\varphi_0)$ the intrinsic bias.

Importantly, we refer to $S_T(\vartheta_0)$ and $B_T(\varphi_0)$ as “exact” biases as we evaluate the sample sums inside these expressions for a fixed $T$. Following the literature, we also derive simplified expressions by replacing the sample sums, suitably scaled, by their asymptotic counterparts. The resulting bias expressions, denoted by $\mathcal{S}_T(\vartheta_0)$ and $\mathcal{B}_T(\varphi_0)$, are referred to as “approximate” and will make the representations in Theorem (ref) and Theorem (ref) more intelligible. The simplification of the exact intrinsic bias $B_T(\varphi_0)$ is based on finding $\mathcal{B}(\varphi_0) = \lim_{T \rightarrow \infty} T B_T(\varphi_0)$ and defining $\mathcal{B}_T(\varphi_0) = \mathcal{B}(\varphi_0)/T$. The simplified expression of the exact score bias $S_T(\vartheta_0)$ is more intricate because its order of magnitude depends on $d_0$. From Theorem (ref), $T S_T (\vartheta_0) = O(1)$ when $d_0 > 1/2$, such that the approximate bias can also be based on finding $\mathcal{S}(\vartheta_0) = \lim_{T \rightarrow \infty} T S_T(\vartheta_0)$ and writing $\mathcal{S}_T(\vartheta_0) = \mathcal{S}(\vartheta_0)/T$. However, when $d_0 < 1/2$, it is shown that $T S_T(\vartheta_0) = O(\log(T))$. Yet, we derive in (ref) the asymptotic expansion $T S_T(\vartheta_0) =\mathscr{S}_T(\vartheta_0) + o(1)$ that allows us to define $\mathcal{S}_T(\vartheta_0) = \mathscr{S}_T(\vartheta_0)/T$. Note that the exact and the approximate biases are potentially very different from each other. For instance, the number of sample sums in $B_T(\varphi_0)$ that are approximated by their limiting value in $\mathcal{B}_T(\varphi_0)$ equals $3 (p+1)^3 + (p+1)^2$. If $T$ is relatively small, the performance of $\mathcal{B}_T(\varphi_0)$ can therefore deteriorate substantially relative to that of $B_T(\varphi_0)$. lieberman2005expansions present a similar argument in the context of Edgeworth expansions of the memory parameter.

It follows from Lemmata (ref), (ref), (ref) and (ref) that analogues of (ref) also hold for $\hat{\vartheta}_{\mu_0}$ and $\hat{\vartheta}_m$. In particular, analogous derivations result in expressions of the score bias $S_T (\vartheta_0)$ involving $ E ( \mathit{D_{\vartheta}L}_{\mu_0}^*(\vartheta_0) )$ and $ E (\mathit{D_{\vartheta}L}_{m}^*(\vartheta_0))$, respectively. Clearly, both expectations are equal to zero, the former due to Theorem (ref) and the latter by construction. The following theorem is the main result of this paper and presents the approximate biases of $\hat{\vartheta}$, $\hat{\vartheta}_{\mu_0}$ and $\hat{\vartheta}_{m}$. The proof of the theorem is given in Appendix (ref) and an extension that discards Assumption (ref) in Theorem (ref) of Section (ref).

The theorem generates five important insights. First, the intrinsic bias can be interpreted as the bias that remains even if the true value of the deterministic component $\mu_0$ were known. Secondly, the intrinsic bias of all three estimators is the same. Thirdly, the memory parameter $d$ solely affects the bias through the score bias, not through the intrinsic bias, the latter merely depending on the short-run dynamics $\varphi_0$. Fourthly, the score bias of the CSS estimator $\hat \vartheta$ dominates its intrinsic bias only in the stationary region because the strength of the constant term is sufficiently strong to impact the stochastic component, see the discussion following Theorem (ref). Fifthly, as johansen2016role note, the key factor to assess the distortion in testing or constructing confidence intervals for $\vartheta_0$ is the relative bias, i.e.\ the ratio of asymptotic bias to asymptotic standard deviation. Taking account of the asymptotic distribution of the three estimators, Theorem (ref) implies that their relative bias is of order $O(T^{-1/2})$ in the non-stationary region. In the stationary region, the relative bias is of order $O(T^{-1/2} \log(T))$ for the CSS estimator with unknown $\mu_0$, while it is of order $O(T^{-1/2})$ for the CSS estimator with known $\mu_0$ and for the MCSS estimator. Thus, especially in the stationary region, a $t$-test for the memory parameter or for the short-run dynamics is more accurate when using the MCSS estimator compared to the CSS estimator. Later in the empirical study, we will exploit this feature to our advantage.

Given the bias expression of $\hat \vartheta_m$ in (ref) the MCSS estimator can be further improved by removing its intrinsic bias. However, $\mathcal{B}_T (\varphi_0)$ depends on the unknown true short-run dynamics $\varphi_0$, implying that a two-step procedure is required for implementation: In the first step, estimate the model using the MCSS estimator to obtain estimates of the short-run dynamics $\hat{\varphi}_m$. Then, in the second step, use these estimates to correct for the intrinsic bias. We call this refined estimator the bias-corrected MCSS $(bcm)$ estimator:

The following corollary shows the validity of this two-step procedure. The proof is given in Appendix (ref).

Note that the correction can of course also be based on the estimated exact intrinsic bias $B_T (\hat \varphi_m)$. That is indeed the version we include in the simulation study in Section (ref).

Special cases

In this section, we apply the general bias expressions of Theorem (ref) to some special cases of interest. In particular, Section (ref) covers the ARFIMA(0,$d$,0) model and compares the results with those of lieberman2005expansions and johansen2016role. Section (ref) focuses on the ARFIMA(1,$d$,0) model, comparing it with nielsen2005finite. Lastly, Section (ref) looks at short-memory models, concluding with the biases of the AR(1) model and a comparison with tanaka1984asymptotic.

ARFIMA(0,\texorpdfstring{\(d\)}{d},0) model

johansen2016role consider the purely fractional model for $d_0 > 1/2$ when the process involves a finite-number $N_0$ of unobserved pre-sample values. They derive the asymptotic biases for $\hat{d}$ and $\hat{d}_{\mu_0}$. Their results for $N_0 = 0$ are effectively included in Theorem (ref) when $\omega(L;\varphi) = 1$. For this special setting, the present section adds the asymptotic biases of $\hat{d}$ and $\hat{d}_{\mu_0}$ when $d_0 < 1/2$, as well as the asymptotic bias of the MCSS estimator $\hat{d}_m$. Theorem (ref) provides a summary. Note that no use is made of Assumptions (ref) and (ref). Note also that $\omega(L;\varphi) = 1$ implies that $\vartheta = d$.

The expressions for the score and intrinsic biases in Theorem (ref) have simplified considerably. In particular, the intrinsic bias, derived by johansen2016role for $d_0>1/2$, is independent of $d_0$ and therefore identical in both regions. The same bias term appears in lieberman2004expansions for the bias of the estimated memory parameter based on the maximum likelihood estimator in the type-I fractional model with $0 < d_0 < 1/2$ and $\sigma^2$ as well as $\mu$ known, see johansen2016role for a discussion. Since $\mathcal{B}_T$ is pivotal, it can now be easily eliminated. The bias-corrected MCSS ($bcm$) estimator defined in Section (ref) is now simply $\hat{d}_{bcm} = \hat{d}_{m} + T^{-1} 3\zeta_3 \zeta_2^{-2} $, with $bias(\hat{d}_{bcm}) = o(T^{-1})$, as can be directly seen from Corollary (ref).

Analysing the theoretical bias terms through numerical comparisons may assist in building up an intuition. Table (ref) therefore presents the theoretical biases, up to $o(T^{-1})$ terms, of the CSS estimator with unknown and known $\mu_0$ and of the MCSS estimator, for selected values of $d_0$ and $T$. It is evident that the bias of the CSS estimator decreases in both the stationary and non-stationary region as $d_0$ decreases, and decreases everywhere as $T$ increases. As johansen2016role note, the bias of $\hat d$ is equal to that of $\hat d_{\mu_0}$ when $d_0 = 1$, since in this case, $\mathcal{S}_T(d_0) = 0$. Yet it is curious to see that bias($\hat{d}$) is actually smaller than bias($\hat{d}_{\mu_0}$) for $d_0 = 1.1$ and $d_0 = 1.2$. In fact, this occurs for all $d_0 > 1$. The reason is that $\mathcal{S}_T(d_0)$ becomes positive for $d_0 > 1$, thereby increasing the negative intrinsic bias $ \mathcal{B}_T$. The proof is omitted, but the result follows from the fact that $\Psi(d_0) - \Psi(2d_0-1)$ is monotonically decreasing for $d_0 > 1/2$ and from abramowitz1964handbook, stating that $\Psi(d_0) = \log(d_0) + O(d_0^{-1})$ for $d_0 \rightarrow \infty$.

As mentioned in Section (ref), johansen2016role generalise the way the process is initialised by allowing for $N_0$ unobserved pre-sample values and by setting $x_t = 0$ only when $t < 1-N_0$, with $N_0$ a fixed non-negative integer. They show that in the purely fractional model with $d_0 > 1/2$, these unobserved pre-sample values introduce an additional bias term for the CSS estimator. We prove in Theorem (ref) in Appendix (ref) that the MCSS estimator $\hat d_m $ also eliminates the bias associated with the unknown-level parameter $\mu$ in this extended model. It does not, however, address the misspecification bias due to the unobserved pre-sample values.

ARFIMA(1,\texorpdfstring{\(d\)}{d},0) model

In a Monte Carlo simulation performed by nielsen2005finite, the ARFIMA(1,$d$,0) model undergoes thorough analysis. This section provide an explanation of their findings. In the following theorem, we present the approximate bias expression for this particular model. An extension that does not make use of Assumption (ref) is contained in Theorem (ref) in Appendix (ref).

Figure (ref)(a) plots the approximate intrinsic bias $\mathcal{B}_T (\varphi_0)$ of both the memory parameter and the autoregressive coefficient as given in (ref). Interestingly, the biases for these two parameters exhibit near-symmetry. Specifically, the bias for the memory parameter tends to be negative, while the bias for the autoregressive parameter is typically positive. Furthermore, both biases are relatively small and almost linear when $\varphi_0$ is below 0. However, beyond this value, the absolute value of the biases grows until reaching a peak at around $\varphi_0 = 0.5$. For larger values of $\varphi_0$, the biases decrease rapidly in absolute value start. This same behaviour was noticed in a Monte Carlo simulation by nielsen2005finite. They noted that the memory parameter's downward bias is particularly pronounced when the AR coefficient is either 0 or 0.4 and that the estimation methods seem robust against stronger positive AR coefficient, such as 0.8. This aligns with the bias expression in Theorem (ref).

Figure (ref)(b) plots the exact intrinsic bias $B_T (\varphi_0)$ of both the memory parameter and the autoregressive coefficient. $B_T (\varphi_0)$ is based on (ref) in Appendix (ref), simplified appropriately. The patterns of the biases closely resemble those of the approximated intrinsic bias in Figure (ref)(a). However, the specific values differ significantly, particularly in the range between 0 and 0.6. This discrepancy suggests that the asymptotic approximation can lead to notable distortion, especially within this range. Consequently, we suggest to use the exact intrinsic bias when correcting for it. As we will observe later, these biases also align with the findings in our simulation in Section (ref): It is suboptimal to use the approximate intrinsic bias when the sample size is small.

Figure (ref) also illustrates the approximate score bias $\mathcal{S}_T (\vartheta_0)$ (see (ref)-(ref)) as well as the exact score bias $S_T (\vartheta_0)$ (based on (ref)) of the ARFIMA(1,$d_0$,0) model with the memory parameter $d_0$ taking on values of $-0.2$, 0.4, and 1. Several observations can be made: First, the score bias tends to be more pronounced in the stationary region as compared to the non-stationary region, which aligns with what we anticipate on the basis of Theorem (ref). Secondly, there exists a noticeable symmetry between the memory parameter and the autoregressive coefficient. Thirdly, a close match is observed, on the whole, between the exact and approximate score biases for $d_0 = -0.2$ and $d_0$ = 1. However, this correspondence breaks down when $d_0 = 0.4$ and the approximate biases are much larger in absolute value than the exact biases. This distortion arises due to the presence of the term $(1-2d_0)^{-1}$ in (ref), which diverges as $d_0$ approaches 0.5. It is, therefore, recommendable to employ the exact score biases in practical empirical applications. robinson2015efficient who, in a panel setting, correct for the score bias of the CSS estimator, also observe that the exact score bias is preferable. Importantly, this bias is inherently eliminated by the MCSS estimator, obviating the need for additional correction.

Short-memory models

Theorem (ref) also covers bias expressions in cases where long memory is absent in the model. The resulting expressions are presented in Theorem (ref) below. They are straightforward to derive, so details of the proof are omitted. It suffices to mention that the matrix $\mathcal{B}_T(\varphi_0)$ is truncated by removing the components related to long memory and by setting $d_0 = 0$ in $\mathcal{S}_T(d_0,\varphi_0)$. It is important to emphasise that our model is not restricted to ARMA models alone; rather, it encompasses a broader category of short-memory models, with ARMA models being just one particular instance. Indeed, any representation of the short-run dynamics that conforms to (ref) is admissable, including models like Bloomfield's exponential model.

An extension of the theorem that does not make use of Assumption (ref) is contained in Theorem (ref) in Appendix (ref).

It is instructive to compare these expressions to the corresponding results for stationary and invertible ARFIMA models in Theorem (ref). In that context, the bias of the short-run dynamics in the stationary region for the CSS estimator is of order $O(T^{-1} \log(T))$. This, however, is not true for the bias of the short-run dynamics of the CSS estimator for stationary and invertible ARMA models which is of order $O(T^{-1})$, as shown above in Theorem (ref). An extension of a stationary ARMA model to an ARFIMA model increases the bias of the short-run dynamics to be of the same order of magnitude as that of the memory parameter. Nevertheless, the MCSS estimator effectively eliminates this bias by construction, leading to a reduction in the bias order of $\hat{\varphi}$ to $O(T^{-1})$ for general ARFIMA or ARMA models.

Based on Theorem (ref), we can deduce the analytic bias of an AR(1) model as an illustration. The following corollary presents the expressions describing the analytic biases of the three estimators. The proof is straightforward and omitted. An generalisation that does not make use of Assumption (ref) is contained in Corollary (ref) in Appendix (ref).

$\mathcal{S}_T(\varphi_0)$ and $\mathcal{B}_T (\varphi_0)$ correspond to the bias expressions in the AR(1) model with and without constant term, respectively, derived by tanaka1984asymptotic. Note that Tanaka {\it (i)} analyses a well-specified model for an AR(1) process that begins in the infinite past and {\it (ii)} uses exact maximum likelihood to estimate the model, thereby making use of all available information. We get the same results because, without unobserved pre-sample values, our model is also well-specified and the MCSS estimator equally makes best use of the available information.

The role of unobserved pre-sample values

Recall that our model in (ref)-(ref) has so far been supplemented by Assumption (ref), i.e.\ by the simplification that not only the observations $x_t$ but also the unobserved error terms $u_t$ and $\varepsilon_t$ are equal to zero for $t \leq 0$. Assumption (ref) reflects the conditional nature of the CSS estimator. Moreover, the assumption ensures that the exposition of the results in Theorems (ref) and (ref) to (ref) is not cluttered by bias terms that are due to unobserved pre-sample values.

We now extend that setup and discuss generalisations of our results to different initialisation schemes. It will become plain that our results carry over to these extended settings, with the sole difference that the unobserved pre-sample values lead to what we call model misspecification bias. As was pointed out before, the modification factor of the MCSS estimator is constructed to eliminate the unknown-level bias. It does not, however, account for misspecification bias. This is indeed what we showed in Theorem (ref) in the context of the ARFIMA(0,$d$,0) model initialised by $N_0$ unobserved pre-sample values, see the remark at the end of Section (ref) above. A similar result is derived in the present section for the ARFIMA($p_1,d,p_2$) model in (ref)-(ref), initialised by an error process $u_t$ that starts in the infinite past. Put differently, Assumption (ref) is now not imposed. This is the setup used by, for instance, hualde2011gaussian and hualde2020truncated. The following theorem extends our main results in Theorem (ref), its proof is contained in Appendix (ref).

The components $\mathcal{B}_T(\varphi_0)$ and $\mathcal{S}_T(\vartheta_0)$ are identical to the corresponding quantities in Theorem (ref). As before, $\mathcal{B}_T(\varphi_0)$ is the intrinsic bias. $\mathcal{S}_T(\vartheta_0)$ is again the score bias due to estimating the unknown level, yet it is now only one part of the overall score bias: All three bias expressions in Theorem (ref) now contain a second score bias term, viz.\ one that is due to the CSS estimator setting all pre-sample values equal to zero. Its form depends on whether or not a constant term is estimated and is denoted by $S_T^{\mathrm{init},\mu}(\vartheta_0)$ or $S_T^{\mathrm{init},\mu_0}(\vartheta_0)$, respectively\footnote{Note that, in our terminology, we call these misspecification biases exact. The approximate misspecification biases (their asymptotic counterparts) are given in Lemma (ref).}. The expression in (ref) is part of the bias of the CSS estimator with unknown level and of the MCSS estimator. Specifically, the term involving $c_t(\vartheta)$ shows how initial values feed into the estimation of $\mu$ and thus add bias. The expression in (ref), on the other hand, applies when the level is known and the contribution to the bias stems only from the way initial conditions enter the squared residuals.

As is shown by hualde2020truncated, unobserved pre-sample values do not alter the first-order asymptotic properties of the CSS estimator, yet it is plain from Theorem (ref) that they do matter for the bias expressions. Our MCSS estimator removes the unknown-level score bias $\mathcal{S}_T(\vartheta_0)$. johansen2016role suggest a method to reduce the misspecification score bias due to unobserved initial values: In their model, the conditional maximum likelihood estimator is biased due to the existence of $N_0$ unobserved pre-sample values; yet the authors show that that bias decreases as the sample size $N$ gets larger, see their Corollary 1. An application of the results by johansen2016role to our MCSS estimator for the case of $N=0$ is also contained in Appendix (ref).

Simulation

In this section we conduct a simulation study of the finite sample properties of the CSS estimators with unknown and known $\mu_0$ in (ref) and (ref), respectively, and of the MCSS estimator in (ref), together with the refined bias-corrected version\footnote{Note that we use the exact intrinsic bias $B_T$ as a correction. Simulation results using the approximate intrinsic bias are available upon request.} thereof in (ref). In particular, we take as our DGP the model in (ref) with $u_t$ an AR(1) model, i.e.\ $u_t = \varphi_0 u_{t-1} + \epsilon_t$ with $\epsilon_t \sim \textit{NID}(0,1)$. We set $u_0 = 0$ in accordance with Assumption (ref). Without loss of generality, we let $\mu_0 = 0$, since all estimators are invariant to the value of $\mu_0$. In all settings covered by our experiment, we generate $x_t$ for $T = 32, 64, 128, 256$. We let the long memory parameter $d_0$ vary. In particular, we set $d_0 = -0.2,-0.1,\ldots, 1.1,1.2$. For the autoregressive parameter we use $\varphi_0 \in \{-0.5,0,0.5 \}$. We compute the estimates using the optimising intervals $d \in [d_0-5,d_0+5]$ and $\varphi \in [-0.9999,0.9999]$. All results are based on 10,000 replications\footnote{All computations in this paper are done using MATLAB 2019a, see MATLAB. The convergence criteria used for numerical optimisation are the default ones. The code for replicating the results in this paper is available on request. }. We use the fractional difference algorithm of jensen2014fast to generate our models, as well as to compute the objective function of the estimators.

Table (ref) and (ref) present the Monte Carlo biases (multiplied by 100) of the memory parameter and the autoregressive parameter, respectively. We also report the percentage increase of the bias of the CSS estimator relative to the bias of the MCSS estimator by $\Delta \% |\text{bias}|$ in the last column of each $T$. In addition, Table (ref) and (ref) present the Monte Carlo MSE (multiplied by 100) of the memory parameter and the autoregressive parameter, respectively. We now summarise the main findings.

We first discuss the cases where $\varphi_0 = -0.5$ and $\varphi_0 = 0$. The addition of an autoregressive component to the purely fractional model considerably increases the bias of the CSS estimator of $d$, especially in the stationary region, cf.\ the theoretical biases of the purely fractional model in Table (ref). The CSS estimator $\hat{d}$ clearly underestimates the true $d_0$, while the $\hat{\varphi}$ overestimates the true $\varphi_0$. The MCSS estimator, however, reduces a large part of this bias. The largest reduction occurs in the stationary region, which is also expected from Theorem (ref). Importantly, the bias of the MCSS estimator and the bias of the CSS estimator with known $\mu_0$ are close to each other, confirming our theoretical findings that the leading bias of $\hat{\vartheta}_m$ is the same as that of $\hat{\vartheta}_{\mu_0}$. The bias-corrected MCSS estimator further improves upon the MCSS and CSS estimators, as is expected from Corollary (ref). Turning to the MSE, the MCSS estimator consistently outperforms the CSS estimator across nearly all cases. The largest improvement occurs again in the stationary region, which is where the largest bias reduction takes place. The CSS estimator with known $\mu_0$ performs the best and outperforms the MCSS estimator, while the biases are somewhat similar, the variance of the former estimator is significantly lower because $\mu_0$ is known. This estimator also outperforms the bias-corrected MCSS estimator. Although the bias of the bias-corrected MCSS estimator is lower, its variance appears to be relatively larger.

Additional comments can be made regarding the case $\varphi_0 = 0$. This scenario is both interesting and realistic. Typically, the autoregressive lags of a regression model are unknown, and a lag selection procedure, such as the method proposed by box1990time or the use of an information criterion as in huang2022consistent, is used to estimate it. It is not unlikely that the lag selection procedure or the information criterion overestimates the number of lags. The case $\varphi_0 = 0$ resembles a situation in which the lag order in the model is larger than that of the DGP. Then, according to our simulation, the CSS estimator is strongly biased when the level parameter term is not known, even in comparison to the case with $\varphi_0 = -0.5$. The MCSS estimator, as opposed to the CSS estimator, significantly reduces the bias. Therefore, lag selection procedures when nuisance parameter are included in the model could be based on the MCSS estimator rather than the CSS estimator, e.g.\ see lee2015model.

Finally, we now discuss the situation where $\varphi_0 = 0.5$. The comments made above are also true for the non-stationary region. In particular, the bias of the MCSS estimator is close to that of the CSS estimator with known $\mu_0$. In the stationary region, however, the two estimators behave differently. baillie2024combining show that the correlation between the MLE estimates of the short-run dynamics and memory parameter is $-0.94$ for this setting, which might explain the difficulty for our estimators in distinguishing between the autoregressive and memory component and the possible theoretical mismatch of our estimators. Nevertheless, the differences become small for $T = 256$. Furthermore, the CSS estimator of performs the worst in terms of the bias, while the other three estimators significantly improve on this estimator. In terms of the MSE, we see that the MSE of the MCSS estimatorof $d$ and that of the CSS estimator with known $\mu_0$ is significantly lower than that of the CSS estimator. As opposed to that, the MSE of the CSS estimator for $\varphi$ is lower than that of the MCSS estimator and also the CSS estimator with known $\mu_0$ when $T = 32$.

In order to better understand the differences in the bias and MSE, we have plotted in Figure (ref) densities of the three estimators of $d$ (left panels) and $\varphi$ (right panels) for $T = 32$ (upper panels) and $T = 256$ (lower panels), when $d_0 = -0.2$ and $\varphi_0 = 0.5$. It can be seen that the CSS estimators strongly underestimate the true $d_0 = -0.2$ and strongly overestimate the true $\varphi_0 = 0.5$ for $T = 32$. This strong bias in the CSS estimator contrasts with less variation. The CSS estimator's poor performance extends somewhat to the case $T = 256$. The MCSS and CSS estimators are well-centred, but this centring comes at the cost of an increase in the variance. This explains the differences in the MSE of the CSS estimator relative to that of the MCSS estimator and the CSS estimator with known $\mu_0$. Also, it appears that for small $T$ the MCSS estimator recentres the memory parameter relatively more than the autoregressive component, explaining the differences with the bias of the CSS estimator with known $\mu_0$. Nevertheless, these plots show that MCSS density estimates are more similar to CSS density estimates with known $\mu_0$ than to CSS density estimates with unknown $\mu$. The good finite sample performance of the MCSS estimator is again evident.

Empirical examples

As an illustration of the results derived in Section (ref), we now present three empirical applications, reconsidering the long-memory modelling of classical datasets. What all three applications have in common is that the datasets consist of short time series of 79 to 171 observations each, warranting the use of our MCSS estimator to correct the small-sample bias of the received estimators.

Post-second World War real GNP

sowell1992modeling conducted a well-known empirical analysis of the long-memory behaviour of U.S.\ post-Second World War quarterly, seasonally adjusted, log real GNP. The data\footnote{We use the data provided by potter1995nonlinear in the JAE Data Archive who mentions Citibase as his source, as does sowell1992modeling. The dataset can be downloaded from \url{https://journaldata.zbw.eu/dataset/a-nonlinear-approach-to-us-gnp}.} comprise observations from 1947:2 to 1989:4 and are displayed in panel (a) of Figure (ref).

sowell1992modeling estimates an ARFIMA(3,$d$,2) type-I model of mean-adjusted first differences using full maximum likelihood (ML), basing the lag order on the Akaike information criterion. He obtains an estimated memory parameter of $-0.59$. However, smith1997fractional assert that Sowell's results are substantially biased and especially the memory parameter is strongly underestimated. They propose a simulation-based bias correction of the profile maximum likelihood (BC-PML) estimator, resulting in $\hat{d} = -0.46$. The BC-PML estimator relies on the assumption that the bias is a linear function in the parameters. However, lieberman2005expansions\footnote{lieberman2005expansions consider the profile plug-in maximum likelihood estimator instead of the profile maximum likelihood estimator for tractability reasons.} show this not to be the case for a simple ARFIMA(0,$d$,0) type-I model. We circumvent this problem by using our MCSS estimator, which does not require the bias to be linear in the parameters. Table (ref) presents the CSS and MCSS estimates of $d$ for the ARFIMA(3,$d$,2) type-II model in (ref), along with the ML estimate of sowell1992modeling and the profile maximum likelihood (PML) as well as BC-PML estimate of smith1997fractional. It can be noted, first, that the CSS estimate is of similar order of magnitude as the maximum likelihood estimates, compare e.g.\ CSS and (P)ML. Secondly, the bias-correction increases both the CSS and PML estimates substantially, cf.\ MCSS and BC-PML. In fact, the CSS estimate is increased by a larger margin than the PML estimate. Thirdly, the type-II estimates are less significant than the type-I estimates, and the significance is reduced by the bias-correction. In conclusion, our results indicate that the long memory parameter is closer to zero than previously thought, even relative to its standard error.

Extended Nelson-Plosser dataset

There is a long-standing controversy on whether it is apt to describe the 14 time series in the well-known nelson1982trends dataset, as extended by schotman1991bayesian\footnote{The dataset can be downloaded from \url{http://korora.econ.yale.edu/phillips/data/np&enp.dat} and is included in the R package `tseries'.}, by unit root processes. More recently, the literature on long memory processes has broadened the debate by considering a fractional integration parameter $d$ that can take any value on the real line instead of merely zero or one. Yet the test statistics for the null hypothesis of $d=1$ tend to be close to their critical values, impeding strong conclusions. Prominent papers are, amongst others, crato1994fractional, gil1997testing, shimotsu2010exact and la2019saddlepoint.

Our enquiry proceeds in two stages: First, we revisit crato1994fractional\footnote{Unfortunately, we did not succeed in replicating the results of crato1994fractional. They use Sowell's Fortran program GQSTRFRAC, which is not available to us. Also, hassler2019time mentions an error in the autocovariance formula of sowell1992modeling. We hence exercise caution in interpreting their results.} who use profile maximum likelihood (PML) to estimate an ARFIMA type-I model. We compare their PML to the CSS and MCSS estimates of $d$ in our type-II setting. Secondly, we conduct unit root tests and relate them to the results obtained in the frequency-domain setting considered by gil1997testing and shimotsu2010exact. This comparison is of interest because the MCSS estimator shares one interesting characteristic with frequency-domain estimators, namely that the leading bias of the estimator is not altered by an inclusion of a level parameter, a feature not present in PML or CSS.

The extended Nelson-Plosser dataset consist of 14 annual macroeconomic series, starting between 1860 and 1909 and running to 1988, and are displayed in Figure (ref). For the analysis, all of the series are log-differenced\footnote{The “differencing and adding back” technique, a commonly used method to simplify estimation by removing drift through differencing, has been found to deliver inconsistent CSS estimates in type-II models when the data in levels exhibit a memory parameter of less than 0. As a solution to this problem, hualde2020truncated recommend modelling the data in levels instead of first-differences or, alternatively, employing a single dummy variable to capture the initial observation. Implementing this latter approach, our results remain qualitatively the same.}, except for the bond yield, which is merely in differences. Table (ref) displays the PML, CSS and MCSS estimates of the memory parameter as well as their respective standard error and the $t$-statistics for testing the unit root null that $d = 0$. Following crato1994fractional, the model selection is based on the BIC\footnote{huang2022consistent have recently shown the BIC criterion to provide consistent selection of the short-run dynamics when based on the CSS estimator in ARFIMA models without constant term.} of the PML estimator. The table reveals that (a) four of the PML $t$-statistics are larger than the 5% critical values of a two-sided test, with a further five being borderline cases, (b) the MCSS estimates are consistently larger than the PML and CSS ones, and (c) of the MCSS $t$-statistics, only three lie in or close to the critical region. Another interesting point to note in Table (ref) is that, for the GNP price deflator, PML and CSS provide a long memory estimate of $-0.39$ and $-0.40$, respectively, while the MCSS estimator yields a value of $+0.22$. This disparity may be attributed to the fact that CSS strongly underestimates the memory parameter when positive AR(1) dynamics are present whereas the MCSS estimator eliminates the bias, as demonstrated in Theorem (ref) and the simulation study presented in Section (ref). In summary, we find greater evidence than in the previous literature in favour of the unit root hypothesis in 11 out of the 14 Nelson-Plosser series.

Let us now turn to the second issue of interest, i.e.\ the comparison of our time-domain estimation to the frequency-domain approaches in gil1997testing and shimotsu2010exact. For the unit root null hypothesis, gil1997testing employ robinson1994efficient's (robinson1994efficient) LM-type test based on the Whittle (W) estimator, while shimotsu2010exact uses a $t$-type statistic based on the extended local Whittle (ELW) objective function. Table (ref) compares the results of the unit root test based on the frequency domain estimators W and ELW with those based on the time domain estimators in Table (ref), a checkmark indicating that $H_0$ is (almost) rejected at the 5% level. Two important observations can be made from this table. First, the tests of gil1997testing and shimotsu2010exact give completely different outcomes, confirming the impression that there is presently no consensus in the literature on the unit root issue. A discussion of the relative merits of the W and ELW estimators is provided in, for instance, hualde2011gaussian. Secondly, the test decisions of gil1997testing are consistent with the majority of the MCSS tests. They only differ for real GNP, the unemployment rate and CPI, the reason for which could be that gil1997testing capture the short-run dynamics solely through AR($k$) components, which may be somewhat restrictive considering that the BIC also discovers MA lags.

Nile data

We now present an empirical application to the classical dataset\footnote{ The dataset used in this analysis can be obtained from the R package `datasets'.} on the annual water flow volume of the Nile for the years 1871 to 1970. The 100 time-series observations are displayed in panel (b) of Figure (ref). Several studies have analysed this dataset either in a long memory or short memory framework, with or without the presence of a break in the time series. hosking1984modeling and boes1989parameter focus on long memory without considering a break. macneill1991search, wu2007inference, macneill2020multiple examine breaks in a short memory time series context. atkinson1997detecting look at breaks in a unit root model. shao2011simple and bet17 address the testing and estimation of a break using a procedure that is robust to long memory although, after identifying a break, they do not proceed to estimating the fractional parameter. In summary, while there appears to be a consensus on including a break in the model, there is disagreement on whether the dynamics are better described by short or long memory. In particular, the literature currently does not consider the estimation of the memory parameter that is robust to a break. This is what we aim to achieve.

To that end, we proceed in two steps: First, we extend our model in (ref) to incorporate a break, i.e.\ $\mu$ in (ref) is replaced by $ \mu_t(\tau) = \mu + \beta I(t \leq \lfloor \tau T \rfloor), $ where the break fraction $\tau \in (0,1)$ is assumed unknown. $\mu_t (\tau)$ can be consistently estimated in a type-II fractionally integrated model with $|d_0| < 1/2$, as shown by chang2016inference and iacone2019testing. It is, however, necessary to generalise our Assumption (ref) such that $q > 1/(1+2d_0)$ moments exist, see johansen2012necessary. In a second step, we employ the filtered observations $\hat{x}_t = x_t - \hat{\mu}_t(\hat{\tau})$ to obtain the CSS estimates $\hat{\vartheta}$ in (ref) and the MCSS estimate $\hat{\vartheta}_m$ in (ref). The consistency of $\hat{\vartheta}$ in this model follows from similar arguments as in robinson2015efficient, that of $\hat \vartheta_m$ in this model is easily obtained because of its asymptotic equivalence to the CSS estimator, see (ref) in Lemma (ref). The model selection procedure suggested by hualde2011gaussian is employed, consisting in a preliminary estimator $\tilde{d}$ of $d$ obtained by local Whittle estimation as in robinson1995gaussian before the procedure by box1990time is applied for selecting the short-run dynamics of $\Delta^{\tilde{d}} \left\{\hat{x}_t\right\}$. The lobato1998nonparametric automatic selection rule of the bandwidth $m$ is used.

In the first step, we find that $\hat \tau = 0.27 $, translating into an estimated break in 1898. This is similar to what most of the aforementioned papers find, and it coincides with the beginning of the construction of the Lower Aswan Dam in 1899. The estimates of the level and break magnitudes of, respectively, $\hat \mu (\hat \tau) = 849.97$ and $\hat \beta (\hat \tau) = 247.78$ imply that the flow volume was reduced by 22%. Note that hosking1984modeling implements an alternative adjustment based on the recommendation of todini1979hydrological, namely that the pre-1903 flows are reduced by 8%.

In the second step, we find a bandwidth of $m = 22$, resulting in a preliminary estimate of $\tilde{d} = -0.05$. The Box-Jenkins procedure indicates that the short-run dynamics are best described by a MA(1) model. The resulting CSS and MCSS estimates are reported in Table (ref), along with their standard errors and $t$-statistics. The results are unambiguous: the MCSS estimate does not provide evidence of long memory in the Nile data once the break is incorporated, with the point estimate of the memory parameter being $-0.12$. The CSS estimate supports this conclusion, with an estimate of $-0.18$. In terms of short-run dynamics, however, CSS and MCSS differ: CSS estimates the MA coefficient to be $0.30$, an estimate that is statistically significant at the 5% level. On the other hand, the MCSS estimate of $0.26$ is insignificant. Given the superior finite sample properties of MCSS, our conclusion is that after incorporating the break, the Nile data is characterised by IID shocks. This finding aligns with that of atkinson1997detecting, supporting their argument that the series can be adequately described by a white noise process once the break is taken into consideration\footnote{ atkinson1997detecting also identifies an outlier in the year 1913. However, even after removing this outlier, our results remain robust.}.

To corroborate our conclusion regarding the memory parameter, we employ the semi-parametric $t$-type statistic in iacone2022semiparametric to test the null hypothesis $H_0 \colon d_0 = 0$ against the alternative hypothesis $H_1 \colon d_0 \neq 0$. This test is designed to be robust against breaks and has the advantage that a parametric specification of the shocks is not needed. The test result, omitted to conserve space, is conclusive and supports our finding: after taking into account the break, there is no evidence that the Nile data exhibits long memory.

Discussion and outlook

Practitioners like the CSS estimator due to its simplicity and effectiveness in estimating both stationary and non-stationary ARFIMA models. Recent work by hualde2020truncated,hualde2021truncated provides the asymptotic justification for using the CSS estimator to estimate models that include deterministic components. However, incorporating a level parameter introduces an additional bias component to the CSS estimator. This bias is due to a biased score which is particularly pronounced when the data are stationary. To address this issue, we propose modifying the CSS profile objective function to create an unbiased score, resulting in a new estimator which we call the modified CSS (MCSS) estimator. This new estimator is straightforward to compute and implement, enabling practitioners to obtain more accurate estimates and less distorted tests and confidence intervals. We illustrate the MCSS estimator by a Monte Carlo simulation and by three classical empirical applications.

Our analysis is for the general ARFIMA($p_1$,$d$,$p_2$) model that includes a constant term and unobserved pre-sample values. Various extensions are conceivable and are of potential interest, yet are beyond the scope of this paper: First, further deterministic components could be included in the model, e.g.\ a linear time trend: Denoting by $X$ a $T$ $\times$ 2 matrix of a constant and a linear trend it can be shown that the modification term for the MCSS objective function turns out to be

This modification term is again simple to calculate. Notably, it is equivalent to that in (ref) if $X$ is only a vector of ones and the degrees of freedom in the exponent are replaced by $T-1$. We expect that analogous modification terms result for more general deterministic components in $X$. We conjecture that the corresponding MCSS estimators improve on the CSS estimators and that their biases are the same as those of the CSS estimators with known deterministic components.

Secondly, while our paper focuses solely on univariate fractional time series, the topic takes on added interest when extended to a panel setting. For instance, robinson2015efficient extend the model presented in equations (ref)-(ref) to a panel framework. While the CSS estimator of $\vartheta$ in a panel setting is consistent under large-$T$ asymptotics, its finite sample properties are deficient due to the presence of fixed effects. To address this issue, the authors propose a bias correction that depends on the true parameters, necessitating the use of estimates to render this correction feasible. However, as the finite sample properties of the CSS estimator is unsatisfactory, replacing the true values by estimates leads to similarly suboptimal bias corrections. As an alternative to improving the small-sample properties of the CSS estimator a similar modification to the CSS objective can be made as in Section (ref). The advantage of this approach is highlighted in the recent work of schumann2023role.

Thirdly, the modification term constructed in this paper is designed to address the unknown-level score bias of the CSS estimator. In principle, our approach can be applied to the elimination of the score bias arising from other sources. Investigating, for instance, whether the misspecification score bias resulting from unobserved pre-sample values can be eliminated by a modification term derived from our principles is a promising direction for future research.

Fourthly, the boundary case of $d_0 = 1/2$ is not covered by our theory. Yet our simulations provide evidence that the behaviour of the MCSS estimator for this parameter value is similar to its behaviour for adjacent parameter ranges. Examining this boundary case analytically is beyond the scope of this paper and therefore left to future work.

\printbibliography