On Least Squares Estimation in a Simple Linear Regression Model with Periodically Correlated Errors: A Cautionary Note

Abstract: In this research the simple linear regression (SLR) model with autocorrelated errors is considered. Traditionally, correlated errors are assumed to follow the autoregressive model of order one (AR(1)). Beside this model we will also study the SLR model with errors following the periodic autoregressive model of order one (PAR(1)). The later model is useful for modeling periodically autocorrelated errors. In particular, it is expected to be useful when the data are seasonal. We investigate the properties of the least squares estimators of the parameters of the simple regression model when the errors are autocorrelated and for various models. In particular, the relative efficiency of those estimates are obtained and compared for the white noise, AR(1) and PAR(1) models. Also, the generalized least squares estimates for the SLR with PAR(1) errors are derived. The relative efficiency of the intercept and slope estimates based on both methods is investigated via Monte-Carlo simulation. An application on real data set is also provided. It should be emphasized that once there are sufficient evidences that errors are autocorrelated then the type of this autocorrelation should be uncovered. Then estimates of model’s parameters should be obtained accordingly, using some method like the generalized least squares but not the ordinary least squares.


Introduction
Regression analysis is a very important statistical method that investigates the relationship between a response variable Y and a set of other variables named as independent variables or predictors X 1 , . . ., X p .An important objective of the built model is the prediction of Y for given values of the predictors.
The simple linear regression (SLR) model is the simplest regression model in which we have only one predictor X.This model, which is common in practice, is written as where Y t , X t are the values of the response and predictor variables in the tth trial, respectively, β 0 and β 1 are unknown parameters and ϵ t are usually assumed to be iid from N(0, σ 2 ϵ ) specially for inference purposes.The terms {ϵ t } satisfying such conditions are named in the time series context as the white noise (WN) process (Wei, 2006).The variable X is usually assumed fixed and non-random.For several predictors, the SLR model generalizes to what is known as the multiple linear regression model.For the estimation of the SLR model there are two common methods, the first is the ordinary least squares (OLS) method, which relies on minimizing the sum of square of errors ∑ ϵ 2 t .For the SLR model (1) the OLS estimators of β 0 and β 1 are It is known that these estimators are unbiased and best linear unbiased estimators (BLUE).
The second method is the maximum likelihood method, which under the assumptions of independence and normality of ϵ t produces again β0 and β1 above (Kutner, Nachtsheim, Neter, and Li, 2005, p. 31-32).
In turn, the fitted SLR model is written as Ŷt = β0 + β1 X t so that the estimated errors or residuals, denoted by e t , are defined as e t = Y t − Ŷt , t = 1, . . ., n. Once, the regression model is fitted, an important step in model building and diagnosis is to check for the assumptions of the model, namely; independence, normality and constant variance of errors.The residuals of the fitted model play a primary role for this purpose.Several graphs and tests for residuals may be used to examine those assumptions.These techniques are usually known as residual analysis.In particular, to examine the independence of error terms we use the plot of residuals against time or order (usually named as the residual plot).Besides, an important test specifically designed for testing the lack of randomness in residuals is the Durbin-Watson test.For a detailed account on various methods for the assessment of the assumptions of regression models, see Kutner et al. (2005).
When the assumptions of the regression model are not valid there exist several remedial measures.In particular, when the error terms are correlated, in which case model ( 1) is named as generalized linear regression (GLR) model, a direct remedial measure is to work with a model that calls for correlated error terms (Kutner et al., 2005, p. 127).A common model of errors in this case is the stationary zero-mean AR(1) model so that ϵ t in (1) now satisfies where u t are assumed iid N(0, σ 2 u ) and |ϕ| < 1.When ϕ = 0, this model reduces to the ordinary SLR model.Other possible models for errors can be selected from the wider class of mixed autoregressive moving average (ARMA) models (Box, Jenkins, and Reinsel, 1994).
In the literature, regression models with autocorrelated errors following some ARMA model have been widely studied (Mohammed and Ibazizen, 2008;Huitema and McKean, 2007;Lee and Lund, 2004;Jeske, Bütefisch, and Song, 1996;Zinde-Walsh and Galbraith, 1991).More general results for regression models with correlated errors are considered, for example, by Grenander (1954), Choudhury, Hubata, and Louis (1999) and Yue and Koreisha (2004).In particular, Ohtani (1990) examined the small-sample properties of the generalized least squares (GLS) estimators and tests of individual regression coefficients with autocorrelated errors.Lee and Lund (2004) considered the OLS and GLS estimation for the linear trend model (3) They have derived explicit formulas for the variances of the OLS and the GLS estimates of β 0 and β 1 with errors following some selected ARMA models.
The GLS estimate of β = (β 0 , β 1 ) ′ under the SLR model ( 1) is where is the design matrix where 1 n is a n-dimensional vector of ones and Γ n is the auto-covariance matrix of the errors ϵ = (ϵ 1 , . . ., ϵ n ) ′ written as where γ k = cov(ϵ t , ϵ t−k ), k = 0, 1, . . . .Besides, the variance of βGLS is given by Grenander (1954) showed that the OLS and GLS estimates are asymptotically equally efficient for SLR with stationary correlated errors.However, once the errors are correlated, the OLS estimates loses its minimum variance property.In a later part of this work, we will see that this fact of deficiency of OLS estimates carries over to the case of periodically correlated errors.
In this article we study the properties of OLS estimates for the parameters of SLR when the errors are periodically correlated.In the time series framework it is found that many real time series exhibit periodic autocorrelations that can not be modelled by ordinary seasonal ARMA models (Tiao and Grupe, 1980;Franses and Paap, 2004;McLeod, 1995).This means that the autocorrelations among successive errors changes from one season to another.For example, assuming that {Y 12k+ν } is a monthly time series then periodic autocorrelation means that corr(Y 12k+ν , Y 12k+ν−1 ) is not constant for ν = 1, . . ., 12. McLeod (1995) suggests a simple graphical method to detect periodic autocorrelations in time series.He also proposed a test of periodically autocorrelated errors.McLeod proposed to apply this test on the residuals resulted from fitting seasonal ARIMA models for some seasonal time series.If significant, this test will indicate that the errors are periodically correlated.

GLR Models with Periodically Correlated Errors
Assume that {Y t } (and possibly {X t }) is a seasonal time series with period ω.Then, if the errors in (1) are correlated then they may inhibit some seasonality.In this case there are several approaches to handle this issue.The first is using seasonal ARMA models (see Box et al., 1994) or one can add some extra regressors to (1) to extract the seasonality as for instance seasonal dummy variables or trigonometric functions.The later approach is useful when the seasonality in {Y t } is deterministic and the former is useful when seasonality is stochastic but homogeneous.
An alternative model that is suitable for modelling seasonality is the periodic ARMA (PARMA) models.Writing the time t in terms of the period ω as kω + ν where ν = 1, . . ., ω denotes the season and k denotes the year, the zero-mean PARMA ω (p(ν), q(ν)) model is written as: ) where {u kω+ν } is a zero-mean WN process with periodic variances σ 2 u (ν), p(ν) is the AR order for season ν and q(ν) is the MA order for season ν and ϕ 1 (ν), . . ., ϕ p(ν) (ν) and θ 1 (ν), . . ., θ q(ν) (ν) are the AR and MA parameters of season ν, respectively (Franses and Paap, 2004).The periodic autoregressive model (PAR) is a special case of the PARMA model.In (7) setting q(ν) = 0 for each ν = 1, . . ., ω we get the equation of the PAR ω (p(ν)) model.For instance, the zero-mean PAR ω (1) model can be written as In fact, this equation can be written as ω equations.For instance, the zero-mean PAR 4 (1) model can be written as PARMA models are not stationary in the ordinary weak sense.They are rather examined for a weaker type of stationarity named as periodic stationarity.This means that the mean and the variance of the time series is constant for each season and periodic with period ω and the autocovariance function depends on the time lag and season only (Ula and Smadi, 1997).For example, the PAR 4 (1 Obeysekera and Salas, 1986).The essence of PARMA models is that it is suitable for modelling periodic correlations.The first motivations for those models started in hydrology then found any applications in economics and other areas (Obeysekera and Salas, 1986;Franses and Paap, 2004).
Therefore, the main idea in this research is to propose the idea of periodically correlated errors in SLR model through using the PAR(1) model, that is with errors following (8).In the next section the properties of the OLS estimates of β 0 and β 1 are to be investigated.Later on, the GLS estimates of β o and β 1 are also derived.

Properties of OLS Estimates with Correlated Errors
The method of moments is one of the most common methods of estimation in statistical inference.It is also common in the context of time series analysis.Assuming that {U kω+ν } is a periodic stationary process, then the seasonal autocorrelation function depends on the time lag and season only and is defined as where γ j (ν) = cov(U kω+ν , U kω+ν−j ) denotes the seasonal autocovariance function (SACVF) and γ 0 (ν) denotes the variance of the process for season ν and time lag j.
Based on an observed realization u 1 , . . ., u mω the moment estimator of ρ j (ν) is where where ūν is the sample mean of data in season ν and m is the number of years of data.
It can be shown that r j (ν) are asymptotically unbiased and consistent estimators of ρ j (ν) (McLeod, 1995).
As far as the PAR ω (1) is considered, it can be proved that the first lag autocorrelations are given by: for ν = 1, . . ., ω.Note that in this case the first order autocorrelation are not the same as AR parameters but a function of them unless all γ 0 (ν) are equal.For a given PAR ω (1) model with known ϕ 1 (ν) and σ 2 u (ν), ν = 1, . . ., ω, then γ 0 (ν) can be obtained, for all ν, by solving the system of equations Then ρ 1 (ν) can be computed using (5).
For the computation of γ j (ν) for PAR ω (1) models the following theorem is useful.
Theorem 1: If {ϵ kω+ν } follows a periodic stationary PAR ω (1) model ( 8), then Proof: It is easy to show that for j = 1, Thus for j = 2, Iterating for k we get (12).When ϵ t is WN we have seen that β0 and β1 are the BLUE and var( β0 ) and var( β1 ) are (Kutner et al., 2005) and Now, the following theorem is useful for the computation of var( β0 ) and var( β1 ) for errors following any stationary stochastic process and in particular the AR(1) model.
Theorem 2: For the GLR model (1) with {ϵ t } being any stationary stochastic process with mean zero and ACVF γ k , then β0 and β1 are unbiased and with Proof: It can easily be shown that where c t and b t are given by ( 17).Thus, Because, ∑ n t=1 c t = 1 and ∑ n t=1 c t X t = 0 then E( β0 ) = β 0 .Similarly we can show that E( β1 ) = β 1 .The results in (15) and ( 16) are direct from the fact that β0 and β1 are linear functions in Y 1 , . . ., Y n with coefficients as in (17).
Theorem 3: For the GLR model (1) with errors following a zero-mean PAR ω (1) model ( 8) with | ∏ 4 ν=1 ϕ 1 (ν)| < 1, then β0 and β1 are unbiased and where Proof: The unbiasedness of β0 and β1 can be proved in a similar manner as Theorem 2. Now we want to prove ( 22) and ( 23).We know that where c t is defined in (17).Now, assuming n = mω where m is the number of years of data, and writing t as (k − 1)ω + ν, then (25) reduces to where I A is the indicator function for the set A. The results for var( β1 ) in ( 23) can similarly be derived.In fact, Theorem 3 is valid for all periodic-stationary processes with seasonal ACVF γ k (ν).

The Relative Efficiency for Correlated Errors
It is known that if θ1 and θ2 are two unbiased estimators of θ, then the relative efficiency of θ1 with respect to θ2 is (Rohatgi, 1984, p. 193) Therefore, for the SLR model as β0 and β1 are unbiased estimators for β 0 and β 1 for errors following WN, AR(1) and PAR(1) models, then we can compare the efficiency of these estimators using (26).For simplicity, the comparison will be carried out with the WN model being the reference case.Those relative efficiencies are easily obtained from the results in the previous section and are explicitly given below.
Corollary 2: The relative efficiency of β0 and β1 for errors following the AR(1) model with respect to the standard SLR model (i.e.WN errors) are ) with K tj and M tj given in ( 20) and ( 21), respectively.
In view of Table 1 we can see that in terms of efficiency, the LS estimators β0 and β1 based on the AR(1) model (b) are more (less) efficient than based on WN model when ϕ < 0 (ϕ > 0).We in fact conclude from this result that when the errors are autocorrelated via an AR(1) scheme, then when the autocorrelation is positive (negative) then the standard formulas for var( β0 ) and var( β1 ) given by ( 13) and ( 14), respectively, do underestimate (overestimate) their true values.Besides, as the autocorrelation becomes stronger the gaps between the estimated and true variances of β0 and β1 increase.Therefore, if the errors are actually autocorrelated but ignoring this fact and using the standard formulas for var( β0 ) and var( β1 ) then all subsequent results, including prediction, are inaccurate (see, for example, Lee and Lund, 2004).When the errors follow the PAR 4 (1) model (c), that is the errors are periodically autocorrelated, then var( β0 ) and var( β1 ) do also differ from their standard formulas in (13) and ( 14).Specifically, it can be seen in Table 2 and 3 that the relative efficiency of β0 and β1 between the WN errors and the PAR 4 (1) model differ from one.However, there is no specific pattern in the relative efficiencies in view of ∏ ϕ 1 (ν) as seen in Table 2. Finally it is worth mentioning that the relative efficiencies for β0 and β1 were very close to each others and have the same direction in all selected cases.Table 2: Relative efficiencies of β0 and β1 for WN/PAR 4 (1).
As the theoretical SACVF of {ϵ t } is unknown, then it is traditionally estimated via the residuals {e t } of the OLS of (1).That is, if C n is the same as Λ n above with C j (ν) Austrian Journal of Statistics, Vol. 41 (2012), No. 3, 211-226 (defined in (10)) in place of γ j (ν) then For the PAR ω (1) error model, γ j (ν) are computed using Theorem 1 where γ 0 (ν) is estimated by C 0 (ν).For the estimation of ϕ 1 (ν), in view of (8), it is obtained through regressing (through origin) of the residuals in the ν-th season on the preceding residuals which belong to the (ν − 1)-th season.
Finally, we notice that although the estimates and standard errors based on OLS and GLS are close to each others, the OLS results here are invalid and those based on GLS should be more accurate.
In Example 2 above the GLS and OLS estimates are compared for some real data.In fact, the implications of the nature and strength of the autocorrelation structure among errors on the relationship between those estimates is a crucial question.From a design point of view, various studies highlight important issues related to effective design under correlation (see for example Müller and Stehlík, 2009).In the following example we study the effect of the autocorrelation design on the estimates of both methods.Albertson, Aylen, and Lim (2002) studied the power of Durbin-Watson test for SLR models with PAR(1) errors.The selected designs in the next example are due to them.
Case 2: ϕ 1 (1) = 0, ϕ 1 (2) = ρ, ϕ 1 (3) = 0, ϕ 1 (4) = ρ, with ρ = −2, −1, −0.5, 0.5, 1, 2. Then for each realization the relative efficiency of the OLS and the GLS estimates of β 0 = 1 and β 1 = 5 is computed.Then the average relative efficiency based on all realizations is obtained and summarized in Tables 4 and 5.In case 1, the geometric mean of the ϕ's is ρ which represents the degree of correlation while τ controls the degree of periodicity.Case 2 represents a semi-periodic PAR(1) model.Notice for case 2 that ρ Table 5: Average relative efficiencies of OLS and GLS estimates for case 2. For case 1, it can be seen that relative efficiency for both intercept and slope are nearly equal.Besides, for negative (positive) autocorrelation, variances due to OLS are larger (smaller) than for GLS.This is true for the semi-PAR(1) model in case 2. In case 1, as the degree of periodicity increases the relative efficiency values, generally increase.

Conclusion and Remarks
In this article we have investigated the properties of the LS estimators of the slope and intercept parameters of the simple linear regression model.We obtained the relative efficiency of those estimators when the errors are correlated via an AR(1) model as well as when the errors are periodically correlated following the PAR(1) model compared to the standard case of WN errors.An example is provided for various models of those types.The main conclusions are that the LS estimators are still unbiased but the ordinary formulas for var( β0 ) and var( β1 ) are not reliable and may underestimate or overestimate the actual variances of those estimators.Those variances are affected by the strength and direction of autocorrelation among errors when the errors follow the AR(1) and PAR(1) models and are also affected by the nature of periodic autocorrelation when the errors follow the PAR(1) model.Therefore, the known deficiency of OLS estimates when the errors are AR(1) also extends to the case of periodically correlated errors.We believe that the conclusions above extend to the multiple regression model.This point can be investigated in a future research.
When the errors are autocorrelated then the ordinary formulas for var( β0 ) and var( β1 ) should not be used.In this case, a suitable alternative is using GLS method which is defined here when the errors are periodically correlated.Some other remedial measures can also be applied such as including some seasonal component in the regression model as dummy variables or trigonometric functions.

Figure 1 :
Figure 1: The time series plot of quarterly U.S. airline passenger miles(in millions).
below.This graph shows a nearly linear trend with seasonality.The fitted OLS of the linear trend model is given by Ŷt = 104.608+ 0.866t, so that

Table 4 :
Average relative efficiencies of OLS and GLS estimates for case 1.