Psicothema was founded in Asturias (northern Spain) in 1989, and is published jointly by the Psychology Faculty of the University of Oviedo and the Psychological Association of the Principality of Asturias (Colegio Oficial de Psicólogos del Principado de Asturias).
We currently publish four issues per year, which accounts for some 100 articles annually. We admit work from both the basic and applied research fields, and from all areas of Psychology, all manuscripts being anonymously reviewed prior to publication.

  • Director: Laura E. Gómez Sánchez
  • Frequency:
         February | May | August | November
  • ISSN: 0214-9915
  • Digital Edition:: 1886-144X
  • Address: Ildelfonso Sánchez del Río, 4, 1º B
    33001 Oviedo (Spain)
  • Phone: 985 285 778
  • Fax: 985 281 374

Psicothema, 2002. Vol. Vol. 14 (nº 3). 669-672


Jaume Arnau y Roser Bono

University of Barcelona

The statistical analysis of short time series designs is influenced by the presence of serial dependence. Therefore, it is important to correctly estimate the first-order autocorrelation in behavioral data. The empirical bias is an indicator generally used to evaluate the adequacy degree of an estimator. This paper presents the Bias program, a Monte Carlo simulation program that generates first-order autoregressive processes and calculates the bias in three autocorrelation estimators (r1, r1+ and r1’) for different values of the lag-one autocorrelation parameter and sample sizes. The program has been designed with MATLAB programming language and it runs on IBM-PC compatible computers with a 486 or later processor.

Programa para el cálculo del sesgo empírico de estimadores de autocorrelación. El análisis estadístico de los diseños de series temporales cortas está influido por la presencia de dependencia serial. De ahí la importancia de estimar correctamente la autocorrelación de primer orden en datos conductuales. El sesgo empírico es un indicador generalmente usado para evaluar el grado de precisión de un estimador. Este artículo presenta el programa Bias, un programa de simulación Monte Carlo que genera procesos autorregresivos de primer orden y calcula el sesgo de tres estimadores de autocorrelación (r1, r1+ y r1’) para diferentes valores del parámetro de autocorrelación de retardo uno y tamaños muestrales. El programa ha sido diseñado mediante el lenguaje de programación MATLAB y funciona en ordenadores IBM-PC con un procesador 486 o superior.


The study of autocorrelation in time series was considered from the end of the seventies (Bono & Arnau, 1996, 2000; Busk & Marascuilo, 1988; Escudero & Vallejo, 2000; Glass, Willson, & Gottman, 1975; Gorsuch, 1983; Hartman et al., 1980; Huitema, 1985, 1988; Jones, Vaught, & Weinrott, 1977; Rosel, Jara, & Oliver, 1999; Sharpley & Alavosius, 1988; Suen, 1987; Suen & Ary, 1987; Vallejo, 1994, among others). In the nineties, the works were oriented towards the correction of the bias generated in short time-series by the conventional autocorrelation estimator r1. Matyas and Greenwood (1991) verified that the difference between the mean of the conventional estimator r1 and the value of the autocorrelation parameter ρ1 increases with short time series. These authors concluded that the relationship between ρ1 and the mean r1 was nonlinear. The greatest coincidence between calculated mean r1 and ρ1 occurred when ρ1 was around -0.2 or -0.3. According to these values, when ρ1 increased positively, the calculated mean r1 underestimated r1 with a negative bias, while when the ρ1 increased negatively, the calculated mean r1 underestimated ρ1, but with a positive bias. Based on this simulation study, it was deduced that when small sample sizes were used, the r1 values were biased. Huitema and McKean (1991) proposed the modified estimator r1+. This estimator corrects the empirical bias for positive values of the autocorrelation parameter by adding 1/n to r1. Arnau (1999) and Arnau and Bono (2001) proposed the estimator r1’ which consists of the correction of r1 by the absolute value of a fitting model. This model is obtained from the polynomial function of the bias for each sample size. The estimator r1’ reduces the bias for small sample sizes, both for positive as well as negative values of ρ1.

This article describes a computer program using MATLAB, Version 5.2 (1998). The program, called Bias, makes it possible to calculate the empirical bias in r1, r1+ and r1’ by Monte Carlo simulation. This program requires the specification of the parameter ρ1 and the sample size (n). The structure of the program is as follows (see the Appendix): simulation of lag-one autoregressive processes, calculation of the estimators r1 and r1+, polynomial models for different sample sizes, calculation of the estimator r1’, and estimation of the empirical bias in r1, r1+ and r1’. Bias runs under the Windows 3.1 or later operating system on IBM-PC compatibles with a 486 or later processor, and 8 MB of extended memory.

In the first place, the program simulates first-order autoregressive processes. The randn function (Forsythe, Malcolm, & Moler, 1977) is used to generate observations from a lag-one autoregressive process:

Yt= ρ1Yt-1 + et (1)

where Yt is the score on the response measure at time t, ρ1 is the autocorrelation parameter and et is a random normal variate with a mean of zero and a variance of one. Each time series starts with a normal variate (Y0) having zero mean and variance 1/(1-ρ21). For each combination of n and ρ1, 10,000 samples are generated. The first 30 observations of each series are dropped to eliminate any dependency between them (DeCarlo & Tryon, 1993; Huitema & Mckean, 1991, 1994a, 1994b, 1994c). After this, the program calculates the estimators r1 (equation 2) and r1+ (equation 3) for each simulated time series. The conventional estimator r1is defined by the following equation:

where Yt and Yt+1 are the data obtained in the intervals t and t+1

and . The estimator r1 is biased with small sample sizes because the numerator only includes n-1 terms instead of n. Huitema and McKean (1991) proposed the modified estimator r1+ that corrects, to some extent, for the negative bias generated by positive autocorrelations

r1+ = r1+ (1/n) (3)

The empirical bias function of the conventional estimator r1 mildly drifts from the linearity (Arnau, 1999; Arnau & Bono, 2001; DeCarlo & Tryon, 1993; Huitema & McKean, 1991, 1994a, 1994b). Models with significant coefficients of the empirical bias function for n= 6, 10, 20 and 30 were derived with polynomial fitting. Within the MATLAB environment, another program that calculated the empirical bias of the estimator r1 in function of ρ1 and n by Monte Carlo simulation was designed (Arnau, 1999; Arnau & Bono, 2001). The polytool function of the Statistics Toolbox of the MATLAB was used for polynomial fitting. With the significant polynomial coefficients at 5% (Table 1) the following fitting models based on n were designed and were introduced into the Bias program.

Modeln=6= -0.1648 – 0.5643r1 – 0.0916r21 (4)

Modeln=10= -0.0972 – 0.3760r1 – 0.0676r21 (5)

Modeln=20= -0.0482 – 0.2028r1 – 0.0333r21 (6)

Modeln=30= -0.0373 – 0.1360r1 (7)

For sample sizes other than those studied in this paper, it is previously necessary to calculate, by simulation, the empirical biases of the conventional estimator. Once the corresponding empirical biases have been found, the polynomial fitting is carried out with polytool function specifying the ρ1 values, empirical biases previously obtained and the degree of the polynomial fit. In this way, the significant coefficients are obtained. These must be included in the polynomial models for different sample sizes statement indicating the value of n.

Using the equations 4-7, the Bias program calculates the estimator r1’ by adding the corresponding polynomial model in absolute values to r1.

r1’= r1 + Fitting Model (8)

Finally, for each simulated sample, the program calculates the empirical bias in r1, r1+ and r1’, which is the difference between the mean value of the estimator and the value of the parameter ρ1.

The Bias program only calculates r1’ for four sample sizes (n= 6, 10, 20 and 30). For other n values, the polynomial fitting of the bias must be performed first in order to determine the corresponding correction model and to introduce it into the program.

Running the program

The Bias program is a function M-file and the input variables in the MATLAB workspace, autoregressive parameter (rho) and sample size (n), must be introduced into it. As can be seen in the Appendix, the first line of the function M-file specifies the file name (bias), the input variables names (rho and n) and output variables names (r1_bias, r1plus_bias and r1prime_bias). The following expression must be specified in the MATLAB workspace in order to run the program: [r1_bias, r1plus_bias, r1prime_bias]= bias(rho,n). The MATLAB executes the commands in Bias program and the results are shown in the MATLAB workspace.


Interested users who prefer not to type the program can request a file of the program by e-mailing the second author ( It has been designed with the Version 5.2 of MATLAB but it can also work with previous versions.


The research was supported by Grant PS95-0228 from the General Council of Scientific and Technical Investigation, resolution of the Secretary of State of Universities and Investigation of the Ministry of Education and Culture (Spain).


Program Listing

function [r1_bias, r1_plus_bias, r1_prime_bias]=bias(rho,n)

% ****************************************************


% ****************************************************

% Number of simulations

for s=1:10000

% Number of observations


% Variance


% Standard deviation


% Series with mean 0 and fixed standard deviation


% Random series with mean 0 and variance 1


% First value of the series


% Next values of the series

for t=2:obs




% Eliminations of the first 30 observations


% **************************************************


% **************************************************












% **************************************************


% **************************************************


% ****************************************************


% ****************************************************

if n==6

model = -[(-0.1648)+r1.*-0.5643+(r1.^2).*-0.0916];

elseif n==10

model = -[(-0.0972)+r1.*-0.3760+(r1.^2).*-0.0676];

elseif n==20

model = -[(-0.0482)+r1.*-0.2028+(r1.^2).*-0.0333];

elseif n==30

model = -[(-0.0373)+r1.*-0.1360];


% **************************************************


% **************************************************


% **************************************************


% **************************************************




Arnau, J. (1999). Bias reduction in the estimation of autocorrelation in short time series. Metodología de las Ciencias del Comportamiento, 1, 25-37.

Arnau, J. & Bono, R. (2001). Autocorrelation and bias in short time series: an alternative estimator. Quality and Quantity, 35, 365-387.

Bono, R. & Arnau, J. (1996). C statistic power analysis through simulation. Psicothema, 8, 699-708.

Bono, R. & Arnau, J. (2000). Designs of small samples: Analysis by generalized least squares. Psicothema, 12 Suppl. 2, 87-90.

Busk, P.L. & Marascuilo, L.A. (1988). Autocorrelation in single-subject research: A counterargument to the myth of no autocorrelation. Behavioral Assessment 10, 229-242.

DeCarlo, L.T. & Tryon, W.W. (1993). Estimating and testing autocorrelation with small samples: a comparison of the C-statistic to a modified estimator. Behavior Research Theory, 31, 781-788.

Escudero, J.R. & Vallejo, G. (2000). Comparison of three alternative methods for the analysis of interrupted time-series. Psicothema, 12, 480-486.

Forsythe, G.E., Malcolm, M.A. & Moler, C.B. (1977). Computer methods for mathematical computations. Englewood Cliffs, NJ: Prentice Hall.

Glass, G.V., Willson, V.L. & Gottman, J.M. (1975). Design and analysis of time series experiments. Boulder: Colorado Associated University Press.

Gorsuch, R.L. (1983). Three Methods for analyzing limited time-series (N of 1) data. Behavioral Assessment, 5, 141-154.

Hartmann, D.P., Gottman, J.M., Jones, R.R., Gardner, W., Kazdin, A.E. & Vaught, R.S. (1980). Interrupted time-series analysis and its application to behavioral data. Journal of Applied Behavior Analysis 13, 543-559.

Huitema, B.E. (1985). Autocorrelation in applied behavior analysis: A myth. Behavioral Assessment 7, 107-118.

Huitema, B.E. (1988). Autocorrelation: 10 years of confusion. Behavioral Assessment 10, 253-294.

Huitema, B.E. & McKean, J.W. (1991). Autocorrelation estimation and inference with small samples. Psychological Bulletin, 110, 291-304.

Huitema, B.E. & McKean, J.W. (1994a). Two reduced-bias autocorrelation estimators: rF1and rF2. Perceptual and Motor Skills, 78, 323-330.

Huitema, B.E. & McKean, J.W. (1994b). Reduced bias autocorrelation estimation: Three jackknife methods. Educational and Psychological Measurement, 54, 654-665.

Huitema, B.E. & McKean, J.W. (1994c). Test of H0: ρ1= 0 for autocorrelation estimators: rF1 and rF2. Perceptual and Motor Skills, 78, 331-336.

Jones, R.R., Vaught, R.S. & Weinrott, M. (1977). Time-series analysis in operant research. Journal of Applied Behavior Analysis 10, 151-166.

MATLAB. (1998). MATLAB (Version 5.2). Natick, MA: The MathWorks, Inc.

Matyas, T.A. & Greenwood, K.M. (1991). Problems in the estimation of autocorrelation in brief time series and some implications for behavioral data. Behavioral Assessment, 13, 137-157.

Rosel, J., Jara, P. & Oliver (1999). Cointegration in multivariate time series. Psicothema, 11, 409-419.

Sharpley, C.F. & Alavosius, M.P. (1988). Autocorrelation in behavioral data: An alternative perspective. Behavioral Assessment 10, 243-251.

Suen, H.K. (1987). On the epistemology of autocorrelation in applied behavior analysis. Behavioral Assessment 9, 113-124.

Suen, H.K. & Ary, D. (1987). Autocorrelation in applied behavior analysis: Mith or reality?. Behavioral Assessment 9, 125-130.

Vallejo, G. (1994). Evaluation of the effects of intervention in designs of time series in the presence of trends. Psicothema, 6, 503-524.

Impact factor 2021:   JCR (WOS): 4.104 (Q1)   |   SJR (Scopus) : 1.308 (Q1)    |  CiteScore 2020: 5,3