Abstract
This paper proposes an improved likelihood-based method to test for first-order moving average in the disturbances of nonlinear regression models. The proposed method has a third-order distributional accuracy which makes it particularly attractive for inference in small sample sizes models. Compared to the commonly used first-order methods such as likelihood ratio and Wald tests which rely on large samples and asymptotic properties of the maximum likelihood estimation, the proposed method has remarkable accuracy. Monte Carlo simulations are provided to show how the proposed method outperforms the existing ones. Two empirical examples including a power regression model of aggregate consumption and a Gompertz growth model of mobile cellular usage in the US are presented to illustrate the implementation and usefulness of the proposed method in practice.
1. Introduction
Consider the general nonlinear regression model defined by where is a -dimensional vector of regressors and is a -dimensional vector of unknown parameters. The regression function is a known real-valued measurable function.
Nguimkeu and Rekkas [1] proposed a third-order procedure to accurately test for autocorrelation in the disturbances, , of the nonlinear regression model (1) when only few data are available. In this paper we show how the same approach can be used to test for first-order moving average, MA, disturbances in such models. We therefore assume in the present framework that the disturbance term in Model (1) follows a first-order moving average process, MA, where the error terms are independently and normally distributed with mean and variance and the process is stationary and invertible with .
Although regression models with moving average disturbances are surely no less plausible a priori than autoregressive ones, they have been less commonly employed in applied research. The main reason for this lack of interest has been argued by some authors, for example, Nicholls et al. [2], as the computational difficulty in estimating these models. However, with the recent improvements in computer capabilities and numerical algorithms, such concerns have been greatly alleviated. Unlike in an autocorrelated error process where all error terms are correlated, the MA error process allows correlation between only the current period and the next period error terms. This feature of the MA error process makes it very attractive in situations where precisely this type of correlation structure is more appropriate and instructive. Many real data applications can be found in Franses [3], Hamilton [4], and Lütkepohl and Kratzig [5].
The method proposed in this paper is based on a general theory recently developed by Fraser and Reid [6] and possesses high distributional accuracy as the rate of convergence is whereas the available traditional asymptotic methods converge at rate . The proposed method is accordingly referred to as third-order method and its feature makes it more convenient and particularly attractive for small sample size models. The approach proposed in this paper is also a general procedure for which the Chang et al. [7] statistics, designed for testing moving average in pure univariate time series, is a special case. More specifically, Chang et al. [7] assumed the mean to be some constant whereas the current paper considers a general form for the mean function to be , where is a vector of unknown parameters and the function is known. Although the parameter of interest exists only in the variance-covariance matrix of the error terms and may seem unrelated to the mean function at first, estimation of in regression analysis provides estimates that depend on both the estimates and the functional form of the mean function. Hence, inference about crucially relies on the nature of . This paper proposes a general procedure that accommodates various forms of the regression functions including nonlinear, linear, or constant functions.
In Section 2, we discuss the nonlinear regression model with MA errors and explain how the -value function for the correlation parameter of interest, , can be obtained. We also show how the proposed method can be adopted for a more general higher-order moving average process. For a thorough exposition of the background methodology of third-order inference, we refer the reader to the discussions in Fraser and Reid [6] and Fraser et al. [8]. Numerical studies including Monte Carlo simulations as well as real data applications on aggregate consumption and mobile cellular usage in the US are presented in Section 3 to illustrate the superior accuracy of the proposed method over existing ones and its usefulness in practice. Finally, some concluding remarks are given in Section 4.
2. The Model and the Proposed Method
In this section we present the nonlinear regression model and describe a procedure that can be used to obtain highly accurate tail probabilities for testing first-order moving average disturbances in this model. We will follow the notation in Nguimkeu and Rekkas [1] (hereafter denoted by NR) as close as possible.
2.1. Model and Estimation
Using NR notation we can rewrite the nonlinear regression model (1) with error terms defined by (2) in a reduced matrix form. Denote by the vector of observations of the dependent variable , the regression vector evaluated at , and the vector of error terms. Then Model (1) with error terms defined by (2) can be rewritten as with where is the tridiagonal matrix defined by
Denote by the vector of parameters of the model. We can then obtain a matrix formulation for the log-likelihood function defined by where is the determinant of . It can be shown (see, e.g., [9]) that
The likelihood function defined in (6) can also be easily obtained by adapting the procedure of Pagan and Nicholls [10] who derived it for the general MA() errors in the linear regression framework. The elements of the positive definite matrix , the inverse of , can be found in Whittle ([11], page 75). The element () of , denoted , is given by
The maximum likelihood estimator of can be calculated by solving simultaneously the first-order conditions, , using an iterative procedure such as Newton’s method. We denote by the matrix of partial derivatives of , where is the gradient of . The first-order conditions for the MLE are given by where and .
Let denote the hessian matrix of , . The information matrix of our model can be obtained by calculating the second-order derivatives of the log-likelihood function: where is the element of the -vector and .
Note that the above derivation assumes that the MLEs are obtained as interior solutions. But in practice, the observed likelihood function could be monotone. For such problematic cases, an adjustment procedure has been proposed by Galbraith and Zinde-Walsh [12]. It is also known that in addition to the common convergence problem related to nonlinear optimization, another issue specific to moving average models is the so-called “pileup” or boundary effect (see Stock [13] for a detailed survey). Nonconvergence can arise in nonlinear optimization because the objective function is “flat” in the neighbourhood of the maximum or because there are multiple local maxima. This problem is usually circumvented in practice by using prior or theoretical information to choose starting values or ordinary least squares estimates if linear model is a special case. Alternatively, a grid search can provide appropriate initial estimates if no prior or theoretical information is available. (In the examples provided in Section 3, we used theoretical information to find initial values in Section 3.2. and ordinary least squares estimates as starting values in Section 3.3. In both cases, convergence is reached after only few iterations.) The other issue specific to moving average models is the existence of singularities at the exact upper and lower limits of the moving average parameter region which may produce apparent maxima outside of the restricted feasible region. One way to deal with this problem is to slightly narrow the limits of the parameter space. In any case, one would need to examine the profile-likelihood (obtained by concentrating out in the likelihood function) and not simply rely on the maximum likelihood estimate provided by the score equations.
In the rest of the paper, we denote by our parameter of interest so that , where the corresponding nuisance parameter is . With these notations, the observed information matrix is then given by where the last display expresses as the inverse of the asymptotic covariance matrix of .
2.2. Likelihood-Based First-Order Approximations Methods
Based on the maximum likelihood estimation described above, two familiar methods can be used for testing : the Wald test and the signed log-likelihood ratio test.
The Wald test is based on the statistic () given by where is the maximum likelihood estimate of and the term is represented in the estimated asymptotic variance of given in (11).
Since is asymptotically distributed as standard normal, an -level test can be performed by comparing with , the percentile of the standard normal. Although the Wald method is simple, it is not invariant to parameterization.
The signed log-likelihood ratio test is given by the statistic (), which is one of the most common inferential methods that is invariant to parametrization, and is defined by where is the constrained maximum likelihood estimate of for a given value of . These constrained MLEs, , can be derived by maximizing the log-likelihood function given by (6) with respect to and for a fixed value . This leads to the following first-order conditions for the constrained maximum likelihood estimator:
Similar to the overall likelihood case, we can construct the observed constrained information matrix using the following second-order derivatives: where is the element of the -vector .
It can be seen from the above expressions that the mean and the variance parameters are observed-orthogonal so that the observed constrained information matrix can be written in a simple block-diagonal form. (This is computationally more convenient than using a full matrix expression. The author thanks one of the anonymous referees for pointing this out).
Note that if the interest is also in inference on the regression parameters , or on the variance parameter , then the above derivations can as well be readily applicable by setting , or , and finding the constrained maximum likelihood estimators as well as the observed constrained information matrix accordingly. In the empirical example given in Section 3.3, we illustrate this by performing inference on all the parameters of the model, that is, both the regression and the innovation parameters.
It is well known that the statistics given in (12) and (13) are first-order methods; that is, they are asymptotically distributed as standard normal with first-order accuracy as their limiting distribution has an rate of convergence. Tail probabilities for testing a particular value of can be approximated using these statistics with the cumulative standard normal distribution function , that is, or . However, the accuracy of these first-order methods suffers from the typical drawbacks of requiring a large sample size.
2.3. The Proposed Method
Various methods have been proposed to improve the accuracy of first-order methods in recent years. In particular, Barndorff-Nielsen [14, 15] showed than an alternative approximation of the -value function of , , can be derived with third-order accuracy using a saddlepoint method: The statistic is known as the modified signed log-likelihood ratio statistic and the statistic is the signed log-likelihood ratio statistic defined in (13). The statistic is a standardized maximum likelihood departure whose expression depends on the type of information available. Approximation (16) has an rate of convergence and is thus referred to as a third-order approximation. Barndorff-Nielsen [14] defined for a suitably chosen ancillary statistic in special classes of models. But Fraser and Reid [6] formally introduced a general third-order approach to determine an appropriate for (16) in a general model context. Their approach can be used for inference on any scalar parameter in any continuous model setting with standard asymptotic properties without the requirement of an explicit form of ancillary statistic. Moreover, the existence of an ancillary statistic is not required, but only ancillary directions are. Ancillary directions are defined as vectors tangent to an ancillary surface at the data point.
To obtain the ancillary directions for our nonlinear regression model define by (3), we consider the following full dimensional pivotal quantity: where the matrix is such that . Such a matrix can be obtained, for example, by Cholesky decomposition of . An expression for the elements of , which is found in Balestra [16], is given by Clearly, the quantity is a vector of independent standard normal deviates. The ancillary directions are then obtained as follows:
Hence, a reparametrization specific to the data point can be obtained as follows: Since the sample space gradient of the likelihood is given by the components of the new locally defined parameter are written as The parameters , , and have dimensions , , and , respectively, so that the overall parameter is a -dimensional vector. To reexpress our parameter of interest on this new canonical parameter scale, we require and .
We have and we also have By change-of-basis, the original parameter of interest can then be recalibrated in the canonical parameter, , scale by
The modified maximum likelihood departure is then constructed in the space, and an expression for is given by where and are the observed information matrix evaluated at and observed nuisance information matrix evaluated at , respectively, calculated in terms of the new reparameterization. The determinants can be computed as follows:
These calculations can then be used to obtain the modified maximum likelihood departure measure formulated in (26). The statistics and can then be used in the Barndorff-Nielsen [14] expression given by (16) to obtain highly accurate tail probabilities for testing a specific value of the parameter of interest. The -value function, , obtained from (16) can be inverted to build a confidence interval for , defined by
It is important to note that although the discussion in this section is provided for a general nonlinear regression model, the proposed methodology can be adopted to accurately test for MA disturbances in linear regression models as well as in pure univariate time series models. For example, if the regression function is linear, the proposed procedure can be applied by setting and replacing and by and , respectively, in the above formulas. If the model is a univariate time series with a scalar mean , then by setting , where 1 is an -column vector of ones, the proposed procedure is equivalent to the Chang et al. [7] approach.
The extension of the proposed procedure to higher-order moving average processes also follows along similar lines. It however requires to appropriately redefine the covariance matrix, , featured in (4). For example, suppose that the error structure is MA(), ; that is, where , are the unknown moving average parameters for different lags, the ’s are as before, and the roots of the polynomial equation , each have modulus less than one. Then, it can be shown (see, e.g., [17]) that the matrix takes the form where is an matrix with at place and elsewhere and As before, the exact log-likelihood function of the model is given by (6) and improved inference on each of the moving average parameters , , can be performed using the above procedure.
3. Numerical Studies
This section presents Monte Carlo simulation studies to evaluate the empirical performance of the values obtained by the testing procedure as well as two empirical examples using real data. Results of the proposed test are compared to two tests that are commonly used for testing nonlinear regression model disturbances: the Wald test given by (12) and the signed log-likelihood ratio departure (hereafter denoted by LR) derived from the log-likelihood function given in (13). Since the proposed test is based on the Barndorff-Nielsen [14] approximations, it is hereafter denoted by BN.
3.1. Monte Carlo Simulations
The goal of the simulation study is to evaluate the small sample properties of the proposed method and compare them with other methods. The accuracy of the different methods is evaluated by computing confidence intervals for in each case and using statistical criteria such as “central coverage,” that is, the proportion of the true parameter value falling within the confidence interval, “upper error,” that is, the proportion of the true parameter value falling above the upper limit of the confidence interval, “lower error” that is, the proportion of the true parameter value falling below the lower limit of the confidence interval, and “average bias,” that is, the average of the absolute difference between the upper and the lower probability errors and their nominal levels. These criteria are standard in the literature and have been considered, for example, by Fraser et al. [18] and Chang and Wong [19].
The setup of the Monte Carlo simulation is the logistic growth model that is widely used in many applications (see [1, 20–23]). The model is defined by with disturbances generated by the MA process defined by
The errors are assumed to be distributed independently as standard normal. The explanatory variable is a vector of deterministic integers running from to . The parameter values are given as , , and . Samples of sizes and are used along with 5,000 replications each. The moving-average coefficient takes values and confidence intervals of 90% and 95% are computed for each value of the moving average parameter. Note that the value corresponds to the null hypothesis of independence in the disturbances.
Only simulation results for the 95% confidence interval are reported. Results for all other simulations performed are available from the author. Table 1 reports simulation results for a sample size of , while Figure 1 graphically illustrates the results when the sample size is .
Based on a 95% confidence interval, it is expected that each error probability is 2.5% with a total of 5% for both errors, and the average bias is expected to be 0. The superiority of the proposed third-order method (BN) is clear, with central coverage probability as high as 95.05% and average bias as low as of 0.0047. Moreover, while the third-order method (BN) produces upper and lower error probabilities that are relatively symmetric, with error probabilities totaling approximately 5-6%, those produced by the first-order methods (LR and Wald) are heavily skewed, and in the case of the Wald test statistic, the total error probability reaches as high as 18.5%. The asymmetry in the upper and lower tail probabilities for both these methods are persistent and can be evidenced in Figure 1 where upper and lower error probabilities are plotted for various values of used in the simulation. The discrepancy between many values of the first-order methods and the nominal value of 0.025 is unacceptably large. The average bias and coverage probability are aslo provided in Figure 1. It can be seen from this figure that the proposed method gives results very close to the nominal values whereas the first-order methods give results that are less satisfactory. On the whole, simulations show that a significantly improved accuracy can be obtained for testing first-order moving average disturbances in nonlinear regression models using a third-order likelihood method.
3.2. Empirical Example 1
This empirical example uses a nonlinear regression to fit a model of mobile phone usage in the US for the period 1990–2010. The response function commonly used to model the diffusion of emerging technology is the Gompertz growth function (see, e.g., [24–26]). The Gompertz curve is particularly attractive for technological modeling and forecasting because of its S-shape that well illustrates the life-cycle of new technological products over time. The model considered is the Gompertz regression of mobile cellular usage defined by where is the number of mobile cellular subscriptions per 100 people at time , the parameter represents the asymptote toward which the number of subscriptions grows, reflects the relative displacement, and the parameter controls the growth rate (in scaling) of the mobile phone subscriptions. The model is estimated using data from the 2012 World Telecommunication Indicators database and covers the period 1990 through 2010 which corresponds to the uptake and diffusion period of mobile cellular use in the US.
Figure 2(a) depicts the scatter plots of mobile phone subscriptions over time in the US for the period 1990–2010 as well as the nonlinear regression fit. As the S-curve theory of new technologies life-cycle predicts, the empirical shape of mobile cellular usage in the US exhibits three stages: incubation, rapid growth, and maturity. The nonlinear regression fit of the data using the Gompertz regression function appears to be graphically satisfactory.
(a)
(b)
The estimation of the model is performed using numerical optimization that requires reasonable start values for the parameters. Unfortunately there are no good rules for starting values except that they should be as insightful as possible so that they are close enough to the final values. Given the growing uptake of the mobile cellular technology, it is reasonable to assume that the number of subscriptions will reach at least 100 per 100 people in the long run. A starting value for the asymptote of around is therefore reasonable. If we scale time so that corresponds to 1990, then substituting into the model and using the value from the data, and assuming that the error is , we have . Finally, returning to the data, at time (i.e., in 1991), the number of subscriptions per people was . Using this value along with previously determined start values and again setting the error to gives . The solution to the optimization problem is then reached after six iterations and an achieved convergence tolerance of .
Results for the naive nonlinear least squares (NLS) estimation of the model are given in Table 2. The relative growth rate is estimated at and the displacement parameter is estimated at . More interestingly, the asymptotic upper bound of the S-curve is estimated at ; that is, in the long run the number of subscriptions per 100 people would be about 131.8. This implies a relatively high rate of mobile cellular usage compared to the population size. A possible explanation is the absence of Mobile Number Portability (MNP) regulations in the US which could result in subscribers obtaining multiple SIM cards in order to enjoy lower on-net call rates.
Figure 2(b) gives the NLS residual plots and Figure 3 plots the autocorrelation function (ACF) and the partial correlation function (PACF) of the residuals. A graphical analysis of the residuals as well as the ACF and PACF plots suggests the presence of first-order positive serial correlation in the error terms of the model. This intuition can be tested by assuming a first-order moving average structure in the disturbances. To test its significance, assume that the error term follows the MA process defined by
(a)
(b)
Maximum likelihood (ML) estimates of the model can then be computed by incorporating the above MA errors assumption. The results of the ML estimation are presented in the right panel of Table 2. Using the algorithm described in Section 2 the solution is reached in sixteen iterations after which any further iteration does not substantially change the results.
It can be seen that compared to the naive nonlinear least square estimation which does not account for serial correlation the maximum likelihood estimation accounting for MA errors produces different results and smaller standard errors for the estimates. For example, the difference between the NLS estimate for and its ML estimate is 0.0121 and this difference is significant. The moving average coefficient is estimated at with a standard error of . The achieved convergence tolerance is .
Based on the MLEs, the -value function for each of the methods (LR, Wald, and BN) is calculated for hypothesized values of and plotted in Figure 4. They report the results of the test on a more continuous scale so that any significance level can be chosen for comparison. The two horizontal lines in the figure are the data lines for upper and lower 2.5% levels. Figure 4 shows a clearly distinct pattern between the third-order method (BN) and the first-order methods (LR and Wald). As an illustration, Table 3 reports the 90% and the 95% confidence interval for from each of the methods. The confidence intervals obtained in Table 3 show considerable differences that could lead to different inferences about . The first-order methods, LR and Wald, give wider confidence intervals compared to the confidence intervals produced from the third-order methods, BN. Theoretically, the latter is more accurate.
The residuals diagnostics depicted in Figure 5 are also instructive. They show standardized residuals from the ML estimation (a) and the autocorrelation function (ACF) of these residuals (b). All the standardized residuals are within and , and the ACF plots show that autocorrelation of any lag is not significantly different from zero. This further supports the MA model as a plausible error structure for the data.
(a)
(b)
3.3. Empirical Example 2
This empirical example uses a power regression to estimate a general specification of the aggregate consumption function. The example is taken from Greene [27]. The model considered is given by where is personal disposable income and is aggregate personal consumption. The Keynesian consumption function is a restricted version of this model corresponding to . With this restriction the model is linear. However, if is allowed to vary freely, the model is a nonlinear regression. The model is estimated using yearly data from the US economy from 1990 to 2010. As with the preceeding example, estimation requires reasonable start values for the parameters. For the present model, a natural set of values can be obtained because, as we can see, a simple linear regression model is a special case. We can therefore start and at their linear least squares values that would correspond to the special case of and hence use as a starting value for . The numerical iterations are therefore started at the linear least squares estimates of , , and for , and convergence is achieved after 8 iterations at a tolerance of .
Results for the naive linear (OLS) and naive nonlinear least squares (NLS) estimation of the model obtained by assuming , as well as the maximum likelihood estimation obtained by incorporating a first-order moving average for the disturbances as in (35), are given in Table 4.
Except for the intercept the estimated nonlinear relationship gives significant estimates from both the proposed and the traditional approaches (see Table 5). In particular, although the power coefficient, , is numerically close to unity, it is statistically different from one so that the nonlinear specification is a valid alternative with this sample of observations. Moreover, the nonlinear specification has a slightly better explanatory power than the linear specification as the goodness-of-fit is estimated at 0.9998 and 0.9973, respectively. On the other hand, the moving average parameter is estimated at , very far from the boundary, implying that the likelihood function is likely well-behaved.
However, the residuals do not appear to display serial correlation as suggested by Figure 6. In fact, both the first-order asymptotic test statistics and the proposed third-order procedure fail to reject the null of independence of the disturbances against the MA structure, as evident by the confidence intervals obtained in Table 5.
(a)
(b)
Although in this application all the methods lead to the same conclusion about parameter significance, it is evident that confidence intervals from the first-order methods (LR and Wald) are nearly indistinguishable while confidence intervals based on third-order method (BN) significantly differ from these intervals.
4. Concluding Remarks
In this paper, a third-order likelihood procedure to derive highly accurate -values for testing first-order moving average disturbances in nonlinear regression models is illustrated. It is shown how the proposed procedure can also be used to test for MA disturbances in linear regression models and in pure univariate time series models as well. Monte Carlo experiments are used to compare the empirical small sample performance of the proposed third-order method with those of the commonly used first-order methods, particularly, likelihood ratio and Wald departures. Based on the comparison criteria used in the simulations, the results show that the proposed method produces significant improvements over the first-order methods and is superior in terms of central coverage, error probabilities, and average bias. The results suggest that using first-order methods for testing MA disturbances in nonlinear regression models with small samples could be seriously misleading in practice. The implementation of the proposed procedure is illustrated with two empirical examples that use the Gompertz growth model to estimate the diffusion of mobile cellular usage in the US over time and the power regression model to estimate aggregate consumption as a function of disposable income. MATLAB code is available from the author to help with the implementation of this method.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.