A comparative simulation study of AR(1) estimators in short time series

Various estimators of the autoregressive model exist. We compare their performance in estimating the autocorrelation in short time series. In Study 1, under correct model specification, we compare the frequentist r1 estimator, C-statistic, ordinary least squares estimator (OLS) and maximum likelihood estimator (MLE), and a Bayesian method, considering flat (Bf) and symmetrized reference (Bsr) priors. In a completely crossed experimental design we vary lengths of time series (i.e., T = 10, 25, 40, 50 and 100) and autocorrelation (from −0.90 to 0.90 with steps of 0.10). The results show a lowest bias for the Bsr, and a lowest variability for r1. The power in different conditions is highest for Bsr and OLS. For T = 10, the absolute performance of all measurements is poor, as expected. In Study 2, we study robustness of the methods through misspecification by generating the data according to an ARMA(1,1) model, but still analysing the data with an AR(1) model. We use the two methods with the lowest bias for this study, i.e., Bsr and MLE. The bias gets larger when the non-modelled moving average parameter becomes larger. Both the variability and power show dependency on the non-modelled parameter. The differences between the two estimation methods are negligible for all measurements.


Introduction
Time series analysis has been valuable for achieving insight into the nature of longitudinal processes. Especially the autoregressive moving average (ARMA) model (Box and Jenkins 1976) has gained enormous popularity in various research areas. The autoregressive part models the serial dependence between consecutive measurements. The moving average part models the serial dependence between consecutive error terms. The ARMA(p, q) model is given by: h j e tÀj þ e t ; e t $ Nð0; r 2 e Þ; ð1Þ with y t the score at time t ðt ¼ 1; 2; . . .; TÞ, l the population mean, / i the autocorrelation for lag i ði ¼ 1; 2; . . .; pÞ, h j the moving average parameter at lag j ðj ¼ 1; 2; . . .; qÞ and e t the residual. One of the simplest versions of the ARMA(p, q) is the AR(1) model: where, for simplicity, the subscript 1 is omitted from /. Several estimation methods have been proposed to estimate the AR(1) model. These estimation methods include closed form estimation methods, such as the r 1 estimator (Yule 1927;Walker 1931;Box and Jenkins 1976), C-statistic (Young 1941) and Ordinary Least Squares (OLS) estimator, and iterative estimation methods, such as frequentist Maximum Likelihood Estimation (MLE) and Bayesian Markov Chain Monte Carlo (MCMC) estimation. The performance of the closed form estimation methods in terms of efficiency have been examined and compared in some simulation studies (Huitema and McKean 1991;DeCarlo and Tryon 1993;Arnau and Bono 2001;Solanas et al. 2010). Generally, in particular for shorter time series (e.g., length T 50), the closed form estimation methods have been shown to have biased autocorrelation estimates and/or high variability. Because the closed form and iterative estimation methods have not been mutually compared so far, it is unclear which estimation methods perform better in terms of having a low bias and variability under relevant conditions for empirical practice. Further, little is known about the robustness of the specific estimation methods towards misspecification of the model. This knowledge is important to optimize a time series research design, and to select a low-variability, low-bias, and robust method for estimating an AR(1) model in empirical practice.
In this paper, we discuss two studies to assess the relative performance of several estimators of the AR(1) model. We focus on short time series, with a length T between 10 and 100. Even though these lengths are relevant, for example in psychological research, they are not thoroughly studied yet for all estimators we compare. For the autocorrelation we use values between -1 and 1, and hence consider stationary time series. Our first study provides the information needed to make an informed choice between the estimation methods for the AR(1) model. To this end, we selected five popular and/or promising estimation methods. In a simulation study, we compare these methods with regard to bias, standard error, the bias of the standard error, the rejection rate for / ¼ 0, the power for / 6 ¼ 0, and the point and 95 % interval estimates.
Our second study focuses on the issue of robustness. Robustness, as used in this paper, is the resilience to misspecification with regard to the number of parameters. The effects of misspecification of the ARMA(1,1), AR(1) and AR(2) model have been studied for the least squares estimator. For an underspecified model, the parameters become more biased when the unspecified parameters are further from zero (Tanaka and Maekawa 1984). Overspecification of the model gives a larger prediction mean squared error for the estimation of the score at y t (Kunitomo and Yamamoto 1985). To study the robustness with regard to misspecification, we use the two estimation methods that showed the lowest bias in the first study. In the misspecification study we generate the data using an ARMA(1,1) model, but estimate the parameters as if the data was generated using the same AR(1) = ARMA(1,0) model as used in Study 1.
In the next section, we describe the selection process of the estimation methods used in this paper, followed by a short introduction to the estimators. Then, we present the design, performance criteria and the results of the first simulation study, which aims at comparing various estimators when applied to short time series following an AR(1) model. We continue with the design and results of the second simulation study, which aims at exploring the effect of underspecifying a short time series following an ARMA(1,1) model as data following an AR(1) model. We conclude our paper with a discussion of the simulation studies and the implications of the results.

Selection of estimation methods
To start, we performed a literature search towards estimation methods for AR(1) models. Our selection criteria for the papers were as follows: (1) it must discuss one or more simulation studies that compare different estimators of the AR(1) model; (2) it must include conditions with less than 50 time points and a range of values between r 1 and 1 for the autocorrelation /. This literature search revealed the five papers shown in Table 1.
The earliest discussed estimator is the r 1 estimator (Walker 1931), as implemented in the Yule-Walker model (Yule 1927;Box and Jenkins 1976). However, since several studies have found that the bias of r 1 for small samples is large, especially for data with a positive autocorrelation, various alternatives were proposed McKean 1991, 1994;DeCarlo and Tryon 1993;Arnau and Bono 2001;Solanas et al. 2010). A selection of these is given by name in Table 1. Note that most alternatives are based on the original r 1 , as can be deduced from the names using 'r' or 'r 1 ' and a sub-or superscript. In general, the modifications of r 1 showed a smaller bias than r 1 itself, but a larger variability of the estimated autocorrelation McKean 1991, 1994;Arnau and Bono 2001;Solanas et al. 2010), except for the estimators r þ 1 and the C-statistic. In direct comparisons between r þ 1 and the C-statistic, it was shown that the C-statistic had a smaller average bias and a smaller average mean square error, thus a smaller variability, over different values of / than the r þ 1 estimation method.  , 6, 7, 8, 9, 10, 15, 20, 50, 100 Bias (emp), MSE, a, power A comparative simulation study of AR(1) estimators… 3 Apart from the modifications of r 1 , another closed form solution may be used. The ordinary least squares (OLS) estimator is used in many different applications, most notably in regression analysis. Since the autocorrelation may be interpreted as a special kind of regression parameter, OLS can be used to find the autocorrelation. In comparisons, the OLS estimator showed a smaller bias than most derivations from the r 1 estimators (Huitema and McKean 1991;Solanas et al. 2010). However, the OLS estimator also showed a slightly larger mean squared error than most r 1 derivations. These comparisons between estimators reveal a bias-variance tradeoff in the autocorrelation estimator.
Two important methods that are not found in the comparisons listed in Table 1, are the frequentist MLE and Bayesian MCMC estimation. Though simulation studies using MLE have been done, those studies did not include the conditions of our primary interest. For example, the studies that considered different ARMA(p, q)-models (Stoica et al. 1986;Pantula and Fuller 1985;Garcia-Hiernaux et al. 2009 ), had no condition with less than 100 time points (Cox and Llatas 1991) or were aimed at examining other parts of the estimation process, such as deciding on which ARMA(p, q)-model to use (Watson and Nicholls 1992). This was the same for papers using Bayesian MCMC estimation. Examples of this are studies that have no systematic comparison using different estimators (Price 2012), use AR(2) models (West and Wilcox 1996) or use lagged cross-correlation (Zhang and Nesselroade 2007). The MLE and Bayesian MCMC have become often-used methods of analysis in different fields and applications.

Estimation methods
In the next paragraphs we will describe the five different estimation methods used in this paper.

The r 1 estimator in the Yule-Walker method
The Yule-Walker method for ARMA models (Yule 1927;Walker 1931;Box and Jenkins 1976) may be the best known estimation method in time series analysis. It uses the r 1 estimator to estimate the lag 1 autocorrelation: where y t is the observed score at time t, ðt ¼ 1; 2; . . .; TÞ and y is the mean score over the T observations. Asymptotically, the autocorrelation function for this series is biased by Àð1 þ 4/Þ=T (Kendall and Ord 1990). This bias has empirically been shown to be as large as -0.73 for T = 6 and / ¼ 0:90 (DeCarlo and Tryon 1993). This empirical bias is surprisingly close to the asymptotic bias of -0.77. To keep the bias within reasonable limits, Box and Jenkins (1976, pp. 32-33) advise a minimum length of 50 time points for a time series.
The standard error of the/ r 1 is calculated as: wherer 2 y is the estimated variance of y t andr 2 e is the estimated variance of e.
In comparison studies, several other proposals were done to replace the r 1 estimator McKean 1991, 1994;Young 1941). One of these, which outperformed the r 1 estimator and some of the other estimators in several studies, was the C-statistic (Young 1941;DeCarlo and Tryon 1993;Solanas et al. 2010).

C-statistic
The C-statistic (Young 1941) compensates the bias of the r 1 estimator by adding a factor tô / r1 as:/ The/ C is asymptotically unbiased. The/ C has been shown to be a better estimator than / r1 for / for short time series and a positive / (DeCarlo and Tryon 1993; Solanas et al. 2010)). However, the bias still remains quite large (e.g., -0.38 for / ¼ 0:60 and r = 5) and the power remains quite low (e.g., 0:09 for / ¼ 0:60 and r = 5) for short time series (Solanas et al. 2010).
The standard error associated with/ C is: which is obviously only dependent on the number of observations.

Ordinary least squares
The ordinary least squares (OLS) for an AR(1) model is: The asymptotic standard error for/ ols is: The OLS estimation is capable of handling non-stationary data under certain restrictions. This means that it is possible to obtain a non-stationary estimate (i.e., j/ ols j [ 1). To identify possible different behaviours, we distinguish two types of OLS analysis results: OLS-A will refer to the complete results, where OLS-S will refer to the results where the non-stationary results are left out.

Maximum likelihood estimation
The iterative Maximum Likelihood Estimation (MLE) used to estimate the autocorrelation, shares asymptotic properties with the OLS estimation (Lütkepohl 1991, p. 368-370). The MLE method uses a collection of algorithms to find the maximum likelihood for a parameter or model (Durbin and Koopman 2012). In this study, we will compute the MLE A comparative simulation study of AR(1) estimators… 5 with the 'Broyden-Fletcher-Goldfarb-Shanno' algorithm (Byrd et al. 1995). An asymptotic standard error for/ mle may be estimated in the same way as for/ r 1 , using Eq. 3. The asymptotic bias for an AR(1) model with population mean assumed to be zero, is À2/=T. For an AR(1) model with the mean estimated, the asymptotic bias is ðÀ3/ þ 1Þ=T (Tanaka 1984).

Bayesian Markov Chain Monte Carlo
The Bayesian MCMC is the only non-frequentist estimation method considered in this paper. Bayesian analysis uses a prior probability distribution for the parameters, set up before the analysis. This is combined with the observed likelihood, as computed from the observed data, to form the posterior probability of the parameters. This posterior probability can be expressed through Bayes' theorem: pð/jYÞ / ðYj/Þpð/Þ. For the Bayesian analyses we will use MCMC sampling to find the combination of parameter values which gives the highest likelihood.
In these simulation studies we will consider two weak informative Bayesian priors. Since we assume stationarity we restrict ourselves to prior distributions with non-zero probabilities for j/j 1. That is, we consider a flat prior, giving all values of / between -1 and 1 an equal probability: for À 1 / 1: Further, we consider the symmetrized reference prior defined by (Berger and Yang 1994), which is specifically tailored to autoregressive processes. The symmetrized reference prior is given as: This symmetrized reference prior gives a higher probability to higher values of j/j and has a narrower posterior distribution and a smaller mean square error than the flat prior or Jeffrey's prior in the case of AR(1) models (Berger and Yang 1994). We will denote these methods as B f and B sr , respectively.

Research design study 1: comparison of estimators
To compare the various estimators for the autocorrelation (/), we simulate according to an AR(1) model (see Eq. 2). In the generation of the data we vary the length of the time series T and the autocorrelation /. For T we use five different sizes, namely 10, 25, 40, 50 and 100. For /, we use an autocorrelation of -0.90 to 0.90 inclusive, taking steps of 0.10. Earlier studies show that there is a difference between the bias for the negative and positive / for several estimators, including r 1 and the C-statistic (DeCarlo and Tryon 1993; Solanas et al. 2010). This indicates that a thorough test is required to include both positive and negative autocorrelations. Finally, the number of replications must be set. All of the studies in Table 1 have a minimum of 10,000 replications per condition. However, a pilot study showed that the maximum standard deviation of the mean/ over 5000-10,000 replications was 0.0007, when T = 10 and / ¼ 0:7, for all estimators. Therefore we use N = 2000 replications per condition. Considering a fully crossed experimental design, this yields 19 9 5 9 5000 = 475,000 simulated data sets.
Across all conditions, l is set to zero and r 2 e to one, which can be done without loss of generality. This results in a standard normal distribution for y t given /.
Priors We performed a small simulation study to decide on the values for the hyperparameters of the priors in our Bayesian analyses. In the model we use, only the prior distributions for l and r e have such hyperparameters. We used 3 conditions, with / ¼ À0:50; 0 and 0.50, using 1000 replications per condition and 2000 iterations per analysis. We set T = 10, since shorter series provide less data, and will therefore be more strongly influenced by the choice of the prior. For l we used a normal prior with mean and standard deviation as given, and for r e we used a c prior with shape and rate as given in the top part of Table 2.
As can be seen in Table 2, the differences in the estimated parameters are small, especially when taking into account the uncertainty added by the small T. As a result, we based our choice of priors on theoretical grounds. To reduce the influence of the priors, we choose our priors close to the distributions used for the data generation: l $ Nð0; 2Þ and r e $ Cð2; 2Þ: Outcome measures For each data set we obtain different estimators: r 1 , C-statistic, OLS, MLE, B f and B sr . To compare the estimators, we consider the bias of the various estimators of /, their empirical standard error, the bias of the estimated standard error, the rejection rate for / ¼ 0, power for / 6 ¼ 0, and the point and 95 % interval estimates of /. All outcome measures are calculated for each condition and each estimation method.

Variability
To compare the variability of the different estimators over the different conditions, we consider two estimators: the empirical standard error and the bias of the estimated standard error. The empirical standard error shows the variability of the/ across replications. The bias of the estimated standard error shows to what extent the standard error estimated by the estimation method, resembles the empirical standard error.

Empirical standard error: SDð/Þ
The empirical standard error of/ is calculated by: where / is the mean estimated/ over all replications within a condition.
A comparative simulation study of AR(1) estimators… 7 Test 5 Test 6 Test 7 Priors used l

Bias of the estimated standard error
For the frequentist estimators, the estimated standard error SEð/Þ is calculated using Eqs. 3, 4 and 5, and for the Bayesian estimation, the estimated standard error is obtained through MCMC. To estimate the expected value of SEð/Þ for each estimator, we compute the average SEð/Þ over all replications within a condition: SEð/Þ: To assess the bias of the estimated standard error with regard to the observed standard error, we substract the observed standard error, SDð/Þ from the mean estimated standard error, SEð/Þ: Bias of SEð/Þ ¼ SEð/Þ À SDð/Þ:

Rejection rate and power
For each estimation method and condition, we compute the empirical probability (EPr) for rejecting H 0 : / ¼ 0, with a ¼ 0:05. In the condition with / ¼ 0, the EPr indicates the rejection rate or actual a, in all other conditions the EPr equals the actual power. For the r 1 , MLE, OLS-S and C-statistic methods, first a p value is obtained using a t-distribution.
Considering the t-statistic for a correlation coefficient: For the OLS-A method, a t test based on the estimated standard error of/ is applied, since the possibility of/ having a higher value than one in absolute value renders the t-statistic for correlations inapplicable: For the Bayesian estimation methods, we consider the percentage of datasets for which the 95 % credible interval (CrI) does not hold zero. For each condition and method, we then calculate the EPr of rejecting H 0 : / ¼ 0 as: forr 1 ; C À statistic; MLE; OLSÀA; OLS À S : EPr ¼ #ðH 0 isrejectedÞ=N; forB f ; B sr : EPr ¼ #ðCrI does not hold 0Þ=N:

Point and interval estimates for /
To illustrate the joint effects of bias and variability we consider the two estimation methods with the smallest bias, using the point and interval estimates of /. As point estimate we use the mean of/ per condition, for the interval estimation we use the mean 95 percentile of the/ over all replications per condition.

Procedure
For the simulations and analyses we use the program 'R' (R Core Team 2014). The C-statistic was computed directly with the basic functions available. For the Yule-Walker, OLS and MLE methods we use the command 'ar' from the software package 'stats'. The Bayesian analyses are done with the program 'Rstan' (Stan Development Team 2014).

Results study 1
The OLS estimator rendered estimates of / that were higher than one in absolute value, and thus non-stationary, as expected. The highest percentage of non-stationary estimates, 15.1 %, was found for the shortest series, r = 10 and the highest autocorrelation, / ¼ 0:90. For r = 10 and / ¼ À0:90 to / ¼ 0:80, up to 6.8 % of the estimates per condition were non-stationary, with higher percentages associated with higher values of j/j. For T = 25 to 50 and / ¼ 0:50 to 0.90 in absolute value, up to 2.3 % of the estimates were nonstationary. However, the difference in the results was quite small. Thus we will discuss only the OLS-A results for the OLS, which includes all measurements, unless the OLS-S shows a strong deviation from OLS-A.
For the Bayesian analysis, non-convergence is expressed in the potential scale reduction factor,R. The potential scale reduction factor shows the ratio of how much the estimation may change when the number of iterations is doubled, with a perfect 1 indicating that no change is expected (Gelman and Rubin 1992;Stan Development Team 2014). For each estimated parameter /, l and r e , less than 0.39 % of the estimates showed aR above 1.02. Furthermore, a maximum of 2:8 %, found for l as estimated with B f , showed aR above 1.01.

Bias
The bias of the six estimators as a function of / for T = 10, 25, and 50 is presented in Fig. 1. The conditions for T = 40 and are not shown due to their uninformative nature: T = 40 yields results highly similar to T = 50, and T = 100 yields results with hardly any differences between the estimators. As can be seen in Fig. 1, the bias becomes smaller as T increases for all methods, which is to be expected. The relation between the bias and / is roughly linear for all methods, being positive for negative values of / and negative for positive values of /. Further, the bias for positive values of / is larger than the bias for their negative counterparts (i.e. -/). This holds for all values of T and for all methods, except for the C-statistic.
With regard to the ordering of the estimation methods, differences are found between negative and positive values of / and between short time series, T = 10, and longer time series, T C 25. For the shortest time series with T = 10, the differences between the methods with regard to bias are strongly dependent on /. For low, negative values of /, the smallest bias is shown by the OLS, MLE and, to a lesser extent, the r 1 . For positive values of /, the smallest bias is shown by the B sr , followed by the B f . The largest bias for T = 10 is associated with the C-statistic for negative values of /, and the r 1 for positive values of /.
For T C 25 and any /, B sr consistently shows the smallest bias. Just as for the shortest series, the largest bias for T C 25 is associated with the C-statistic for negative values of /, and with the r 1 for positive values of /.

Variability
With regard to variability, the results for T ! 40 are highly similar to the results for T = 25 with regard to pattern of the variability and the order of the estimation methods. The only difference is the decline in absolute size. This prompted us to only explicitly show the results for T = 10 and T = 25 for the empirical standard error and the bias of the estimated standard error.

Empirical standard error: SDð/Þ
The empirical standard error (SDð/Þ) as a function of / is shown in Fig. 2 for T = 10 (panel a) and T = 25 (panel b). For all frequentist estimators, the SDð/Þ for positive values of / is larger than the SDð/Þ for their negative counterparts (i.e. -/), implying that the variability is higher for positive values of / than for negative values of /. For the Bayesian estimators, this differs between values of T and j/j.
With regard to the ordering of the estimation methods for the SDð/Þ, small differences are found between the short time series, T = 10, and longer time series, T ! 25. For the shortest time series, T = 10, and / below 0.40, the lowest SDð/Þ is shown by r 1 , for / above 0.40 this is shown by B f . The highest SDð/Þ for / below 0.30 is shown by B sr , for / above 0.30 this is shown by the OLS estimator. The OLS and MLE stand out due to the continuing increase in the SDð/Þ for higher values of /. For T C 25, the B sr shows an distinct pattern. The B sr shows the lowest SDð/Þ for / below -0.80 and above 0.80, but the highest SDð/Þ for / between -0.6 and 0.60. The lowest SDð/Þ for / between -0.70 and 0.40 is shown by the r 1 . The highest SDð/Þ for / below -0.70 is shown by the C-statistic, for / above 0.70 this is shown by the OLS followed by the MLE. When T increases, the empirical standard error of the different methods become smaller and more similar to each other. For T = 100, the size of the SDð/Þ is between 0.05 and 0.10 for all values of / and all estimators.

Bias of the estimated standard error
The bias of SEð/Þ as a function of / is shown in Fig. 2 for T = 10 (panel c) and for T = 25 (panel d). In general, the bias of SEð/Þ decreases when T becomes larger, indicating a smaller difference between the estimated and the empirical standard errors. For T = 100, the bias of SEð/Þ is between -0.01 and 0.04 for all values of / and all estimators. With regard to the ordering of the estimation methods, small differences are found between T = 10 and longer time series. Differences were also found for different values of /.
The direction of the bias of SEð/Þ differs between the methods and the value of /. For  For T = 10, the smallest bias of SEð/Þ for / below -0.70 is shown by the OLS method. The smallest bias of SEð/Þ for / above -0.50 is shown by the C-statistic, closely followed by the OLS and the MLE. The largest bias of SEð/Þ for / above -0.80, is shown by r 1 , which is joined in this regard by B f and B sr for / between -0.20 to 0.60.
The bias of SEð/Þ for T C 25 is smaller than the bias of SEð/Þ for T = 10 and the different methods are closer together. The domain of / for which the C-statistic shows the largest bias of all estimators increases when T becomes larger; for T = 10 this is when / is below -0.50 and above 0.70, for T = 100 this is when / is below -0.30 and above 0.20. The other estimators show the same pattern and order in the bias of SEð/Þ as for T = 10.

Rejection rate and power
The EPr of the different methods is presented in Fig. 3 for T = 10, 25 and 50, where the EPr at / = 0 indicates the empirical rejection rate and the EPr at / 6 ¼ 0 the empirical power. As with the bias, the EPr for T = 40 and T = 100 are not shown due to their uninformative nature. When looking at the rejection rate, the empirical rejection rate approaches the nominal a as the length of the time series increases, as to be expected. The rejection rate for T = 10 is between 0.01 and 0.04 and for T C 25 between 0.03 and 0.05, for most estimators. The only exception is the rejection rate for OLS-A at T = 40, which is 0.06. At T = 100, the MLE, B sr , OLS-S and B f show a rejection rate of 0.050, which is equal to the nominal a of 0.05. For all practical purposes, the difference in rejection rates between estimation methods is negligable. The power of the estimated / shows a positive relation to the size of T and the absolute value of /, as expected. When we would consider a minimal power of 0.80, for T = 10 this is only found for the estimators OLS and MLE, and at very low values of /, i.e. / À 0:90. For T = 25 and negative /, the power is above 0.80 for / À 0:60 for all estimators except for the C-statistic, which has a power above 0.80 for / À 0:70; for positive values of /, the B sr shows a power above 0.80 for / ! 0:60, for the other estimators this is for / ! 0:70. For larger T, the power reaches 0.80 at lower values of /; for T = 100, the power is 0.80 for j/j ! 0:30.
The order of the estimation methods with regard to the power is consistent for the different values of T. The highest power for negative / is shown by the OLS, for positive / this is B sr . The lowest power for negative / is shown by the C-statistic, for positive / this is /. In general, the difference in power between the methods becomes smaller as T becomes larger.

Point and interval estimates for /
The mean/ with a 95 % estimation interval for the MLE and B sr estimations can be seen in Fig. 4. We only present this for the MLE and B sr , since these methods show the lowest bias. As can be seen in Fig. 4, the 95 % intervals are larger for smaller values of T, indicating a larger variability in the/, as would be expected. The strongest decrease in both variability and bias is from T = 40 to T = 25: for B sr the bias decreases with up to 89 % and the variability with 29-51 %, for the MLE the bias decreases by 82 % and the variability by 34-51 %.

Conclusion
In general, it may be concluded that the Bayesian B sr and the frequentist MLE perform best in terms of bias, but not in terms of variability. With regard to empirical variability, the r 1 performs best. For the bias of the estimated variability, the MLE performs best. Furthermore, the B sr is favorable with regard to power for positive /, showing only slight differences with the OLS estimator for a negative B sr . This leads us to continue with the MLE and the B sr estimators for the misspecification study of this paper.

Research design study 2: robustness
To study the robustness of the estimation methods, we misspecify the model. The data is still analysed as if they stem from an AR(1) model, but we generate the data using an ARMA(1,1) model. We generate data sets for two different sizes of T, namely 25 and 50. For / and h, we use parameters of -0.90 to 0.90 inclusive, taking steps of 0.15. Every condition consists of 5000 replications. Considering a fully crossed design, this yields 13 9 13 9 2 9 5000 = 1,690,000 datasets.
Again, across all conditions, l is set to zero and r 2 e to one. We consider the same outcome measures for Study 2 as we did for Study 1.

Results study 2
We successively present the results on the bias, empirical standard error, bias of the estimated standard error, rejection rate, power, and point and 95 % interval estimates. Note that when h is zero, the simulated data follows an AR(1) model, rendering the results equal to the results discussed in the first study of this paper, apart from small deviations resulting from simulation variability.
As with the first study, we checked theR of the estimated parameters /, l and r e of B sr . For each of the parameters, less than 0:14 % showed anR above 1.02, and less than 1:69 % showed anR above 1.01.

Bias
In Fig. 5, heatmaps for the bias of B sr and MLE for T = 25 and T = 50 are presented, expressing the bias depending on the combination of / and h. The h influences the bias in two ways: first, the bias is smaller when h is close to zero, second, the bias gets larger when h is further from /.
The bias is also influenced by the value of T and the estimation method. When looking at T, in the MLE the bias for T = 50 is larger than the bias for T = 25, unless both h and / are negative. The difference between the bias of T = 50 and the bias of T = 25 ranges for MLE from -0.03 to 0.12 per condition. For the B sr , the bias for T = 50 is smaller than for T = 25, unless / has a value above 0.30. The difference between the bias of T = 50 and T = 25 ranges for B sr from -0.07 to 0.03 per condition. Comparing the estimation methods reveals that the bias is small, and in general slightly larger for the B sr than for the MLE, with the largest difference being 0.11 for / ¼ 0:60, h ¼ 0, and T = 25. The difference between the estimation methods is larger for T = 25 than for T = 50.

Variability
Close inspection of the results for the variability and EPr for the B sr and MLE estimators and the two lengths of T, revealed that the patterns are very similar across methods and different lengths of T. This prompted us to only present the results of B sr and T = 25 in Fig. 6. However, we discuss any quantitative differences between the methods. For comparison purposes, we also plotted the SDð/Þ, the bias of SEð/Þ and the EPr for B sr and T = 25 of Study 1.

Empirical standard error: SDð/Þ
The empirical standard error (SDð/Þ), for B sr with T = 25 and h ¼ À0:45; 0:00, and 0.45, can be seen in Fig. 6 (panel a). Some differences between the SDð/Þ over different values of h, / and T are found. First, the SDð/Þ shows a positive slope over / for negative values of h, and a negative slope over / for positive values of h. Second, the SDð/Þ is smaller for T = 50 compared to T = 25. For the MLE estimator the SDð/Þ is up to 0.07 smaller for T = 50, for the B sr the SDð/Þ is up to 0.08 smaller for T = 50. When comparing methods of estimation, the SDð/Þ of the B sr is larger than the SDð/Þ of the MLE, for both values of T. For T = 25, this differs up to 0.03, for T = 50 this differs up to 0.01 per condition.

Bias of the estimated standard error
The bias of SEð/Þ for B sr with T = 25 and h ¼ À0:45; 0:00, and 0.45, can be seen in Fig. 6  (panel b). The bias of SEð/Þ for most combinations of / and h, where h 6 ¼ 0, is positive and higher than the bias for the AR(1) data. Only for low values of jhj combined with high values of j/j, the bias of SEð/Þ is negative. The bias of SEð/Þ is larger for T = 25 than for T = 50, for all methods and conditions, with differences up to 0.02 for both methods. Furthermore, the bias of SEð/Þ is slightly larger for the B sr estimation than for the MLE, with a maximum absolute difference of 0.01 for T = 25 and 0.03 for T = 50.

Rejection rate and power
The EPr for B sr with T = 25 and h ¼ À0:45; 0:00, and 0.45, can be seen in Fig. 6 (panel c). When h 6 ¼ 0, the curve of the EPr shifts relative to the curve of the AR(1) data. For a negative h, the curve shifts to the right, for a positive h, the curve shifts to the left. In both cases, the amount the curves shifts is roughly equal to the absolute value of h. This is the same for both methods, with the shape of the curve being dependent on T, as in Study 1. The differences between the methods are negligible.

Point and interval estimates for /
The mean/ with a 95 % estimation interval for the MLE and B sr estimations are presented in Fig. 7. As can be seen in Fig. 7, the 95 % intervals are larger for smaller values of T, indicating a larger variability in the/. For both methods, the decrease in the 95 % estimation interval is between 33 and 44 % from T = 25 to T = 50 for the different conditions. In most conditions, the 95 % estimation interval is larger and the mean/ is slightly higher for the B sr than for the MLE.

Conclusion
We found that the further the h deviates from zero, the larger the difference between the / and/ is. The B sr shows a larger bias than the MLE when h is further from zero, showing a larger influence of the h parameter in the estimated /. Furthermore, the observed variability is slightly smaller for the MLE, with the difference between B sr and MLE being larger for T = 25 than for T = 50.

Discussion
We compared five estimation methods for the autocorrelation: the r 1 , C-statistic, ordinary least squares, maximum likelihood estimation and Bayesian MCMC estimation. For the Bayesian MCMC estimation we used both a flat prior and a symmetrized reference prior, giving a total of six autocorrelation estimators. We compared these estimators with regard to bias, variability, rejection rate, power and point and 95 % estimation interval estimates. After this comparison, we selected the Bayesian estimation with symmetrized reference prior and the maximum likelihood estimator to use in a second study. In this small study we assessed the robustness of the methods against underspecification.
The results we found in the first study largely complied with the results from previous studies. For the bias for positive values of /, we found the bias of the C-statistic and the OLS to be smaller than the bias of r 1 , as did previous studies (DeCarlo and Tryon 1993; Huitema and McKean 1991;Solanas et al. 2010). For the empirical standard error, we found smaller values for r 1 than for OLS, as did Huitema and McKean (1991). The low rejection rate we found for the r 1 and the C-statistic confirms to earlier studies (Huitema and McKean 1991;DeCarlo and Tryon 1993;Solanas et al. 2010). The power we found for Fig. 7 Mean of/ (points) with a 95 % percentile interval (lines) for different values of / and T the C-statistic is similar to the power found by Arnau and Bono (2001). However, the results of Solanas et al. (2010) with regard to power were only partly confirmed by our study: for negative / we indeed found a higher power for OLS followed r 1 , than for the C-statistic. But for positive /, we found a higher power for the C-statistic than for r 1 .
The first study showed a strong improvement in all measures for all methods between T = 10 and T = 25. This improvement continued, be it not as strong, for higher values of T. When comparing methods, B sr showed the smallest bias. For the frequentist methods, this was MLE followed by the C-statistic. The smallest empirical standard error is shown by r 1 , the smallest bias of the estimated standard error is shown by the C-statistic, the OLS and the MLE. We further found that a small bias is often paired with a high empirical standard error. The power was rather low for all methods at the lengths of time series we considered. For T = 25, the power is below 80 % for all methods for / between -0.5 and 0.5, for T = 100, the power is below 80 % for / between -0.2 and 0.2. The differences between methods with regard to power are negligible. In research areas where effect sizes are small, this may pose a problem. Some studies use moving windows to assess the stability of parameter estimates over time. For these moving windows, these results indicate that a moving window of at least 50 time points is advisable, especially when the differences in parameters over time are small.
The first study was conducted to explore the differences between estimation methods for the autocorrelation in a single subject design. However, this is only a small step in a large research area. The next step may be to explore these results in multilevel or group analyses, thus when there is not one but multiple subjects per dataset. Another issue may be how the different methods respond to non-stationary data, i.e., j/j [ 1.
In the second study, the robustness of the MLE and B sr to underspecification was examined. In general, we confirmed the notion that the further the unmodelled parameter is from zero, the larger the influence of this parameter is on/ (Tanaka and Maekawa 1984). As with the first study, the empirical standard error decreased when T became larger. However, the bias reacted differently for both methods: for the MLE, the bias became slightly smaller for most conditions, where the bias of B sr became slightly larger for a larger T. The difference in performance for all measurements between the MLE and the B sr was small for both values of T. It was shown that the bias, variability, rejection rate and power were all highly dependent on the value of the non-modeled parameter in the data, h. This can be related to the fact that the autocorrelation of the error also has an influence on the autocorrelation of the total score.
The robustness study was rather small and specific, looking into only one possible way to misspecify the ARMA (1, 1) model. More options within misspecification should be explored to find how robust the estimation methods are with regard to under-, over-and misspecification. Important points are the influence of a misspecified error distribution or overspecification of the model.
In conclusion, we found that the best performing methods for autocorrelation estimation are the Bayesian estimator with symmetrized reference prior and the maximum likelihood estimator. The difference in performance between these two is, for all practical purposes, negligible. The results for the measurements improving greatly between T = 10 and T = 25 and continue to do so, but in a less spectacular fashion. For the misspecification study, we found the size of h, the non-modelled parameter, to be vital for the performance of the estimation methods. The differences between lengths of the series and estimation methods was of lesser influence on the results.