# Improved Inference for Moving Average Disturbances in Nonlinear Regression Models

Authors: Pierre Nguimkeu
Journal: Journal of Probability and Statistics
 (2014)

---

## 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.


---

## Body

## 1. Introduction

Consider the general nonlinear regression model defined by(1)yt=g(xt;β)+ut,
where xt=(xt1,…,xtk)′ is a k-dimensional vector of regressors and β=(β1,…,βk)′ is a k-dimensional vector of unknown parameters. The regression function g(·) is a known real-valued measurable function.Nguimkeu and Rekkas [1] proposed a third-order procedure to accurately test for autocorrelation in the disturbances, ut, 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(1), disturbances in such models. We therefore assume in the present framework that the disturbance term ut in Model (1) follows a first-order moving average process, MA(1),
(2)ut=ϵt+γϵt-1,
where the error terms ϵt are independently and normally distributed with mean 0 and variance σ2 and the process ut is stationary and invertible with |γ|<1.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(1) error process allows correlation between only the current period and the next period error terms. This feature of the MA(1) 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 O(n-3/2) whereas the available traditional asymptotic methods converge at rate O(n-1/2). 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 g(xt,β), where β is a vector of unknown parameters and the function g(·,·) 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 g(xt,β). This paper proposes a general procedure that accommodates various forms of the regression functions g(xt,β) including nonlinear, linear, or constant functions.In Section2, we discuss the nonlinear regression model with MA(1) errors and explain how the p-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 byy=(y1,…,yn)′ the vector of observations of the dependent variable yt, g(x,β)=(g(x1,β),…,g(xn,β))′ the regression vector evaluated at x=(x1,…,xn)′, and u=(u1,…,un)′ the vector of error terms. Then Model (1) with error terms defined by (2) can be rewritten as
(3)y=g(x,β)+u
with
(4)u~N(0,σ2W),
where W is the tridiagonal matrix defined by
(5)W=W(γ)=(1+γ2γ0⋯0γ1+γ2γ⋯00γ1+γ2⋱⋮⋮⋮⋱⋱γ00⋯γ1+γ2).Denote byθ=(γ,β′,σ2)′ the vector of parameters of the model. We can then obtain a matrix formulation for the log-likelihood function defined by
(6)l(θ;y)=-n2log⁡σ2-12log⁡⁡|W|-12σ2{y-g(x,β)}′W-1{y-g(x,β)},
where |W| is the determinant of W. It can be shown (see, e.g., [9]) that
(7)|W|=(1-γ2n+2)(1-γ2).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(q) errors in the linear regression framework. The elements of the positive definite matrix Σ=W-1, the inverse of W, can be found in Whittle ([11], page 75). The (j,i)th element (j,i=1,…,n) of Σ(γ), denoted κji, is given by
(8)κji=(-γ)j-i(1-γ2i)(1-γ2(n-j+1))(1-γ2)(1-γ2(n+1))j≥i.The maximum likelihood estimatorθ^=(γ^,β^′,σ^2)′ of θ can be calculated by solving simultaneously the first-order conditions, lθ(θ^)=0, using an iterative procedure such as Newton’s method. We denote by ∇βg(x,β)=(∇βg(x1,β),…,∇βg(xn,β)) the k×n matrix of partial derivatives of g(x,β), where ∇βg(xt,β)=((∂g(xt,β))/(∂β1),…,(∂g(xt,β))/(∂βk))′ is the k×1 gradient of g(xt;β). The first-order conditions for the MLE θ^ are given by
(9)lγ(θ^;y)=(n+1)γ2n+11-γ2n+2-γ1-γ2-12σ^2×{y-g(x,β^)}′Σ^γ{y-g(x,β^)}=0,lβ(θ^;y)=1σ^2∇βg(x,β^)Σ^{y-g(x,β^)}=0,lσ2(θ^;y)=-n2σ^2+12σ^4{y-g(x,β^)}′Σ^×{y-g(x,β^)}=0,
where Σ^=Σ(γ^) and Σ^γ=(∂Σ(γ)/∂γ)|γ=γ^.Let∇ββ′g(xt,β) denote the k×k hessian matrix of g(xt,β), t=1,…,n. The information matrix jθθ′(θ^) of our model can be obtained by calculating the second-order derivatives of the log-likelihood function:
(10)lγγ(θ^;y)=(n+1)[(2n+1)γ2n+γ4n+2](1-γ2n+2)2-1+γ2(1-γ2)2-12σ^2{y-g(x,β^)}′Σ^γγ{y-g(x,β^)},lγβ(θ^;y)=1σ^2∇βg(x,β^)Σ^γ{y-g(x,β^)},lγσ2(θ^;y)=-12σ^4{y-g(x,β^)}′Σ^γ{y-g(x,β^)},lββ(θ^;y)=-1σ^2∇βg(x,β^)Σ^∇βg(x,β^)+1σ^2∑t=1n‍∇ββ′g(xt,β^)ωt^,lβσ2(θ^;y)=-12σ^4∇βg(x,β^)Σ^{y-g(x,β^)},lσ2σ2(θ^;y)=n2σ^4-1σ^6{y-g(x,β^)}′Σ^{y-g(x,β^)},
where ω^t is the tth element of the n-vector ω^=Σ^{y-g(x,β^)} and Σ^γγ=(∂2Σ(γ)/∂γ2)|γ=γ^.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 σ2 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 θ′=(γ,β′,σ2)=(ψ,λ′), where the corresponding nuisance parameter is λ′=(β′,σ2). With these notations, the observed information matrix jθθ′(θ^) is then given by
(11)jθθ′(θ^)=[-lγγ(θ^)-lγλ′(θ^)-lγλ′(θ^)-lλλ′(θ^)]=[jψψ(θ^)jψλ′(θ^)jψλ′(θ^)jλλ′(θ^)]=[jψψ(θ^)jψλ′(θ^)jψλ′(θ^)jλλ′(θ^)]-1,
where the last display expresses jθθ′(θ^) 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 (q) given by
(12)q(ψ)=(ψ^-ψ){jψψ(θ^)}-1/2,
where θ^=(ψ^,λ^′)′ is the maximum likelihood estimate of θ and the term jψψ(θ^) is represented in the estimated asymptotic variance of θ^ given in (11).Sinceq is asymptotically distributed as standard normal, an α-level test can be performed by comparing |q| with zα/2, the (1-α/2)100th 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 (r), which is one of the most common inferential methods that is invariant to parametrization, and is defined by
(13)r(ψ)=sgn⁡(ψ^-ψ)[2{l(θ^)-l(θ^ψ)}]1/2,
where θ^ψ=(ψ,λ^ψ′)′ is the constrained maximum likelihood estimate of θ for a given value of ψ. These constrained MLEs, θ^ψ=(ψ,λ^ψ′)′=(γ,β^ψ,σ^ψ2), can be derived by maximizing the log-likelihood function given by (6) with respect to β and σ2 for a fixed value ψ=γ. This leads to the following first-order conditions for the constrained maximum likelihood estimator:
(14)lβ(θ^ψ;y)=1σ^ψ2∇βψg(x,β^ψ)Σ{y-g(x,β^ψ)}=0,lσ2(θ^ψ;y)=-n2σ^ψ2+12σ^ψ4{y-g(x,β^ψ)}×Σ{y-g(x,β^ψ)}=0.Similar to the overall likelihood case, we can construct the observed constrained information matrixj^λλ′(θ^ψ) using the following second-order derivatives:
(15)lββ(θ^ψ;y)=-1σ^ψ2∇βg(x,β^ψ)Σ∇βg(x,β^ψ)+1σ^ψ2∑t=1n‍∇ββ′g(xt,β^ψ)·ω~t,lβσ2(θ^ψ;y)=-1σ^ψ4∇βg(x,β^ψ)Σ{y-g(x,β^ψ)},lσ2σ2(θ^ψ;y)=n2σ^ψ4-1σ^ψ6{y-g(x,β^ψ)}Σ{y-g(x,β^ψ)},
where ω~t is the tth element of the n-vector ω~=Σ{yt-g(xt;β^ψ)}.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βs,s=1,…,k, or on the variance parameter σ2, then the above derivations can as well be readily applicable by setting ψ(θ)=βs,s=1,…,k, or ψ(θ)=σ2, 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 O(n-1/2) 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, Φ(q) or Φ(r). 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 p-value function of ψ, p(ψ), can be derived with third-order accuracy using a saddlepoint method:
(16)Φ(r*),wherer*(ψ)=r(ψ)-1r(ψ)log⁡(r(ψ)Q(ψ)).
The statistic r* is known as the modified signed log-likelihood ratio statistic and the statistic r is the signed log-likelihood ratio statistic defined in (13). The statistic Q is a standardized maximum likelihood departure whose expression depends on the type of information available. Approximation (16) has an O(n-3/2) rate of convergence and is thus referred to as a third-order approximation. Barndorff-Nielsen [14] defined Q 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 Q 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 directionsV for our nonlinear regression model define by (3), we consider the following full dimensional pivotal quantity:
(17)z(y;θ)=R{y-g(x,β)}σ,
where the matrix R is such that R′R=Σ. Such a matrix can be obtained, for example, by Cholesky decomposition of Σ. An expression for the elements rij of R, which is found in Balestra [16], is given by
(18)rij=γi-j(1-γ2j)(1-γ2i)(1-γ2(i+1)),j≥i.
Clearly, the quantity z(y;θ) is a vector of independent standard normal deviates. The ancillary directions are then obtained as follows:
(19)V={∂z(y;θ)∂y′}-1{∂z(y;θ)∂θ′}|θ^=(1σ^{y-g(x,β^)}R^-1R^γ{y-g(x;β^)};-∇β′g(x;β^);-1σ^{y-g(x,β^)}).Hence, a reparametrizationφ(θ) specific to the data point y=(y1,…,yn)′ can be obtained as follows:
(20)φ(θ)′=[φ1(θ);φ2(θ);φ3(θ)]=∂l(θ,y)∂y·V.
Since the sample space gradient of the likelihood is given by
(21)∂l(θ,y)∂y=-1σ2{y-g(x,β)}′Σ,
the components of the new locally defined parameter φ(θ)are written as(22)φ1(θ)=-1σ2{y-g(x,β)}′ΣR^-1R^γ{y-g(x;β^)},φ2(θ)=1σ2{y-g(x,β)}′Σ∇β′g(x;β^),φ3(θ)=1σ2σ^2{y-g(x,β)}′Σ{y-g(x,β^)}.
The parameters φ1(θ), φ2(θ), and φ3(θ) have dimensions 1×1, 1×k, and 1×1, respectively, so that the overall parameter φ(θ) is a (k+2)-dimensional vector. To reexpress our parameter of interest on this new canonical parameter scale, we require φθ(θ) and ψθ(θ^ψ).We have(23)φθ(θ)=∂φ(θ)∂θ′=[∂φ1(θ)∂γ∂φ1(θ)∂β′∂φ1(θ)∂σ2∂φ2(θ)∂γ∂φ2(θ)∂β′∂φ2(θ)∂σ2∂φ3(θ)∂γ∂φ3(θ)∂β′∂φ3(θ)∂σ2],
and we also have
(24)ψθ(θ)=∂ψ(θ)∂θ′=[∂ψ(θ)∂γ∂ψ(θ)∂β′∂ψ(θ)∂σ2]=[101×k0].
By change-of-basis, the original parameter of interest can then be recalibrated in the canonical parameter, φ(θ), scale by
(25)χ(θ)=ψθ(θ^)φθ-1(θ^ψ)φ(θ).The modified maximum likelihood departure is then constructed in theφ(θ) space, and an expression for Q is given by
(26)Q=sgn⁡(ψ^-ψ)|χ(θ^)-χ(θ^ψ)|{|j^φφ′(θ^)||j^(λλ′)(θ^ψ)|}1/2,
where j^φφ′ and j^(λλ′) 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:
(27)|j^φφ′(θ^)|=|j^θθ′(θ^)||φθ′(θ^)|-2,|j^(λλ′)(θ^ψ)|=|j^λλ′(θ^ψ)||φλ′(θ^ψ)φλ′(θ^ψ)|-1.These calculations can then be used to obtain the modified maximum likelihood departure measureQ formulated in (26). The statistics Q and r 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 p-value function, p(ψ), obtained from (16) can be inverted to build a (1-α)100% confidence interval for ψ, defined by
(28)[min⁡⁡{p-1(α2),p-1(1-α2)},max⁡⁡{p-1(α2),p-1(1-α2)}].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(1) 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 g(x,β)=x′β and replacing ∇βg(x,β) and ∇ββ′g(xt,β) by x′ and 0, respectively, in the above formulas. If the model is a univariate time series with a scalar mean β, then by setting g(x,β)=1β, where 1 is an n-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,W, featured in (4). For example, suppose that the error structure is MA(q), q>1; that is,
(29)ut=∑s=0q‍γsϵt-s,withγ0=1,
where γs, s=1,…,q are the unknown moving average parameters for different lags, the ϵt’s are as before, and the roots of the polynomial equation ∑s=0q‍γszq-s=0, each have modulus less than one. Then, it can be shown (see, e.g., [17]) that the matrix W takes the form
(30)W=∑s=0q‍ρsLs,
where Ls is an n×n matrix with 1 at (t,t+s) place and 0 elsewhere and
(31)ρs={γs+γs+1γ1+⋯+γqγq-sfors=0,1,…,q,0fors>q.
As before, the exact log-likelihood function of the model is given by (6) and improved inference on each of the moving average parameters γs, s=1,…,q, can be performed using the above procedure.


## 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 byy=(y1,…,yn)′ the vector of observations of the dependent variable yt, g(x,β)=(g(x1,β),…,g(xn,β))′ the regression vector evaluated at x=(x1,…,xn)′, and u=(u1,…,un)′ the vector of error terms. Then Model (1) with error terms defined by (2) can be rewritten as
(3)y=g(x,β)+u
with
(4)u~N(0,σ2W),
where W is the tridiagonal matrix defined by
(5)W=W(γ)=(1+γ2γ0⋯0γ1+γ2γ⋯00γ1+γ2⋱⋮⋮⋮⋱⋱γ00⋯γ1+γ2).Denote byθ=(γ,β′,σ2)′ the vector of parameters of the model. We can then obtain a matrix formulation for the log-likelihood function defined by
(6)l(θ;y)=-n2log⁡σ2-12log⁡⁡|W|-12σ2{y-g(x,β)}′W-1{y-g(x,β)},
where |W| is the determinant of W. It can be shown (see, e.g., [9]) that
(7)|W|=(1-γ2n+2)(1-γ2).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(q) errors in the linear regression framework. The elements of the positive definite matrix Σ=W-1, the inverse of W, can be found in Whittle ([11], page 75). The (j,i)th element (j,i=1,…,n) of Σ(γ), denoted κji, is given by
(8)κji=(-γ)j-i(1-γ2i)(1-γ2(n-j+1))(1-γ2)(1-γ2(n+1))j≥i.The maximum likelihood estimatorθ^=(γ^,β^′,σ^2)′ of θ can be calculated by solving simultaneously the first-order conditions, lθ(θ^)=0, using an iterative procedure such as Newton’s method. We denote by ∇βg(x,β)=(∇βg(x1,β),…,∇βg(xn,β)) the k×n matrix of partial derivatives of g(x,β), where ∇βg(xt,β)=((∂g(xt,β))/(∂β1),…,(∂g(xt,β))/(∂βk))′ is the k×1 gradient of g(xt;β). The first-order conditions for the MLE θ^ are given by
(9)lγ(θ^;y)=(n+1)γ2n+11-γ2n+2-γ1-γ2-12σ^2×{y-g(x,β^)}′Σ^γ{y-g(x,β^)}=0,lβ(θ^;y)=1σ^2∇βg(x,β^)Σ^{y-g(x,β^)}=0,lσ2(θ^;y)=-n2σ^2+12σ^4{y-g(x,β^)}′Σ^×{y-g(x,β^)}=0,
where Σ^=Σ(γ^) and Σ^γ=(∂Σ(γ)/∂γ)|γ=γ^.Let∇ββ′g(xt,β) denote the k×k hessian matrix of g(xt,β), t=1,…,n. The information matrix jθθ′(θ^) of our model can be obtained by calculating the second-order derivatives of the log-likelihood function:
(10)lγγ(θ^;y)=(n+1)[(2n+1)γ2n+γ4n+2](1-γ2n+2)2-1+γ2(1-γ2)2-12σ^2{y-g(x,β^)}′Σ^γγ{y-g(x,β^)},lγβ(θ^;y)=1σ^2∇βg(x,β^)Σ^γ{y-g(x,β^)},lγσ2(θ^;y)=-12σ^4{y-g(x,β^)}′Σ^γ{y-g(x,β^)},lββ(θ^;y)=-1σ^2∇βg(x,β^)Σ^∇βg(x,β^)+1σ^2∑t=1n‍∇ββ′g(xt,β^)ωt^,lβσ2(θ^;y)=-12σ^4∇βg(x,β^)Σ^{y-g(x,β^)},lσ2σ2(θ^;y)=n2σ^4-1σ^6{y-g(x,β^)}′Σ^{y-g(x,β^)},
where ω^t is the tth element of the n-vector ω^=Σ^{y-g(x,β^)} and Σ^γγ=(∂2Σ(γ)/∂γ2)|γ=γ^.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 σ2 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 θ′=(γ,β′,σ2)=(ψ,λ′), where the corresponding nuisance parameter is λ′=(β′,σ2). With these notations, the observed information matrix jθθ′(θ^) is then given by
(11)jθθ′(θ^)=[-lγγ(θ^)-lγλ′(θ^)-lγλ′(θ^)-lλλ′(θ^)]=[jψψ(θ^)jψλ′(θ^)jψλ′(θ^)jλλ′(θ^)]=[jψψ(θ^)jψλ′(θ^)jψλ′(θ^)jλλ′(θ^)]-1,
where the last display expresses jθθ′(θ^) 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 (q) given by
(12)q(ψ)=(ψ^-ψ){jψψ(θ^)}-1/2,
where θ^=(ψ^,λ^′)′ is the maximum likelihood estimate of θ and the term jψψ(θ^) is represented in the estimated asymptotic variance of θ^ given in (11).Sinceq is asymptotically distributed as standard normal, an α-level test can be performed by comparing |q| with zα/2, the (1-α/2)100th 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 (r), which is one of the most common inferential methods that is invariant to parametrization, and is defined by
(13)r(ψ)=sgn⁡(ψ^-ψ)[2{l(θ^)-l(θ^ψ)}]1/2,
where θ^ψ=(ψ,λ^ψ′)′ is the constrained maximum likelihood estimate of θ for a given value of ψ. These constrained MLEs, θ^ψ=(ψ,λ^ψ′)′=(γ,β^ψ,σ^ψ2), can be derived by maximizing the log-likelihood function given by (6) with respect to β and σ2 for a fixed value ψ=γ. This leads to the following first-order conditions for the constrained maximum likelihood estimator:
(14)lβ(θ^ψ;y)=1σ^ψ2∇βψg(x,β^ψ)Σ{y-g(x,β^ψ)}=0,lσ2(θ^ψ;y)=-n2σ^ψ2+12σ^ψ4{y-g(x,β^ψ)}×Σ{y-g(x,β^ψ)}=0.Similar to the overall likelihood case, we can construct the observed constrained information matrixj^λλ′(θ^ψ) using the following second-order derivatives:
(15)lββ(θ^ψ;y)=-1σ^ψ2∇βg(x,β^ψ)Σ∇βg(x,β^ψ)+1σ^ψ2∑t=1n‍∇ββ′g(xt,β^ψ)·ω~t,lβσ2(θ^ψ;y)=-1σ^ψ4∇βg(x,β^ψ)Σ{y-g(x,β^ψ)},lσ2σ2(θ^ψ;y)=n2σ^ψ4-1σ^ψ6{y-g(x,β^ψ)}Σ{y-g(x,β^ψ)},
where ω~t is the tth element of the n-vector ω~=Σ{yt-g(xt;β^ψ)}.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βs,s=1,…,k, or on the variance parameter σ2, then the above derivations can as well be readily applicable by setting ψ(θ)=βs,s=1,…,k, or ψ(θ)=σ2, 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 O(n-1/2) 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, Φ(q) or Φ(r). 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 p-value function of ψ, p(ψ), can be derived with third-order accuracy using a saddlepoint method:
(16)Φ(r*),wherer*(ψ)=r(ψ)-1r(ψ)log⁡(r(ψ)Q(ψ)).
The statistic r* is known as the modified signed log-likelihood ratio statistic and the statistic r is the signed log-likelihood ratio statistic defined in (13). The statistic Q is a standardized maximum likelihood departure whose expression depends on the type of information available. Approximation (16) has an O(n-3/2) rate of convergence and is thus referred to as a third-order approximation. Barndorff-Nielsen [14] defined Q 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 Q 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 directionsV for our nonlinear regression model define by (3), we consider the following full dimensional pivotal quantity:
(17)z(y;θ)=R{y-g(x,β)}σ,
where the matrix R is such that R′R=Σ. Such a matrix can be obtained, for example, by Cholesky decomposition of Σ. An expression for the elements rij of R, which is found in Balestra [16], is given by
(18)rij=γi-j(1-γ2j)(1-γ2i)(1-γ2(i+1)),j≥i.
Clearly, the quantity z(y;θ) is a vector of independent standard normal deviates. The ancillary directions are then obtained as follows:
(19)V={∂z(y;θ)∂y′}-1{∂z(y;θ)∂θ′}|θ^=(1σ^{y-g(x,β^)}R^-1R^γ{y-g(x;β^)};-∇β′g(x;β^);-1σ^{y-g(x,β^)}).Hence, a reparametrizationφ(θ) specific to the data point y=(y1,…,yn)′ can be obtained as follows:
(20)φ(θ)′=[φ1(θ);φ2(θ);φ3(θ)]=∂l(θ,y)∂y·V.
Since the sample space gradient of the likelihood is given by
(21)∂l(θ,y)∂y=-1σ2{y-g(x,β)}′Σ,
the components of the new locally defined parameter φ(θ)are written as(22)φ1(θ)=-1σ2{y-g(x,β)}′ΣR^-1R^γ{y-g(x;β^)},φ2(θ)=1σ2{y-g(x,β)}′Σ∇β′g(x;β^),φ3(θ)=1σ2σ^2{y-g(x,β)}′Σ{y-g(x,β^)}.
The parameters φ1(θ), φ2(θ), and φ3(θ) have dimensions 1×1, 1×k, and 1×1, respectively, so that the overall parameter φ(θ) is a (k+2)-dimensional vector. To reexpress our parameter of interest on this new canonical parameter scale, we require φθ(θ) and ψθ(θ^ψ).We have(23)φθ(θ)=∂φ(θ)∂θ′=[∂φ1(θ)∂γ∂φ1(θ)∂β′∂φ1(θ)∂σ2∂φ2(θ)∂γ∂φ2(θ)∂β′∂φ2(θ)∂σ2∂φ3(θ)∂γ∂φ3(θ)∂β′∂φ3(θ)∂σ2],
and we also have
(24)ψθ(θ)=∂ψ(θ)∂θ′=[∂ψ(θ)∂γ∂ψ(θ)∂β′∂ψ(θ)∂σ2]=[101×k0].
By change-of-basis, the original parameter of interest can then be recalibrated in the canonical parameter, φ(θ), scale by
(25)χ(θ)=ψθ(θ^)φθ-1(θ^ψ)φ(θ).The modified maximum likelihood departure is then constructed in theφ(θ) space, and an expression for Q is given by
(26)Q=sgn⁡(ψ^-ψ)|χ(θ^)-χ(θ^ψ)|{|j^φφ′(θ^)||j^(λλ′)(θ^ψ)|}1/2,
where j^φφ′ and j^(λλ′) 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:
(27)|j^φφ′(θ^)|=|j^θθ′(θ^)||φθ′(θ^)|-2,|j^(λλ′)(θ^ψ)|=|j^λλ′(θ^ψ)||φλ′(θ^ψ)φλ′(θ^ψ)|-1.These calculations can then be used to obtain the modified maximum likelihood departure measureQ formulated in (26). The statistics Q and r 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 p-value function, p(ψ), obtained from (16) can be inverted to build a (1-α)100% confidence interval for ψ, defined by
(28)[min⁡⁡{p-1(α2),p-1(1-α2)},max⁡⁡{p-1(α2),p-1(1-α2)}].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(1) 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 g(x,β)=x′β and replacing ∇βg(x,β) and ∇ββ′g(xt,β) by x′ and 0, respectively, in the above formulas. If the model is a univariate time series with a scalar mean β, then by setting g(x,β)=1β, where 1 is an n-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,W, featured in (4). For example, suppose that the error structure is MA(q), q>1; that is,
(29)ut=∑s=0q‍γsϵt-s,withγ0=1,
where γs, s=1,…,q are the unknown moving average parameters for different lags, the ϵt’s are as before, and the roots of the polynomial equation ∑s=0q‍γszq-s=0, each have modulus less than one. Then, it can be shown (see, e.g., [17]) that the matrix W takes the form
(30)W=∑s=0q‍ρsLs,
where Ls is an n×n matrix with 1 at (t,t+s) place and 0 elsewhere and
(31)ρs={γs+γs+1γ1+⋯+γqγq-sfors=0,1,…,q,0fors>q.
As before, the exact log-likelihood function of the model is given by (6) and improved inference on each of the moving average parameters γs, s=1,…,q, can be performed using the above procedure.


## 3. Numerical Studies

This section presents Monte Carlo simulation studies to evaluate the empirical performance of thep 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
(32)yt=θ1[1+eθ2+θ3xt]-1+ut,t=1,…,n
with disturbances ut generated by the MA(1) process defined by
(33)ut=ϵt+γϵt-1,t=…,-1,0,1,….The errorsϵt are assumed to be distributed independently as standard normal. The explanatory variable xt is a vector of deterministic integers running from 0 to n-1. The parameter values are given as θ1=56, θ2=2.9, and θ3=-0.24. Samples of sizes n=10 and n=20 are used along with 5,000 replications each. The moving-average coefficient takes values γ∈{-0.9;-0.6;-0.3;0;0.3;0.6;0.9} and confidence intervals of 90% and 95% are computed for each value of the moving average parameter. Note that the value γ=0 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. Table1 reports simulation results for a sample size of n=10, while Figure 1 graphically illustrates the results when the sample size is n=20.Table 1
Simulation results forn=10, 95% CI.
γ
Method
Central  coverage
Upper  error
Lower  error
Average  bias
−0.9
LR
0.9120
0.0685
0.0195
0.0245
Wald
0.9155
0.0405
0.0440
0.0173
BN
0.9385
0.0285
0.0330
0.0058
−0.6
LR
0.9285
0.0655
0.0060
0.0298
Wald
0.9250
0.0560
0.0190
0.0185
BN
0.9390
0.0275
0.0365
0.0070
−0.3
LR
0.9180
0.0610
0.0210
0.0200
Wald
0.9015
0.0690
0.0295
0.0243
BN
0.9405
0.0260
0.0335
0.0047
0
LR
0.8510
0.0585
0.0905
0.0495
Wald
0.8150
0.0805
0.1045
0.0675
BN
0.9365
0.0280
0.0355
0.0068
0.3
LR
0.8650
0.0760
0.0590
0.0425
Wald
0.8230
0.0825
0.0945
0.0635
BN
0.9380
0.0295
0.0325
0.0060
0.6
LR
0.9090
0.0370
0.0540
0.0205
Wald
0.8730
0.0410
0.0860
0.0385
BN
0.9360
0.0305
0.0335
0.0070
0.9
LR
0.8890
0.0170
0.0940
0.0385
Wald
0.8630
0.0215
0.1155
0.0470
BN
0.9345
0.0295
0.0370
0.0083Figure 1
Graphical illustration of the simulation results forn=20, 95% CI.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 Figure1 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
(34)yt=β1exp⁡[-β2β3t]+ut,t=0,1,…,n,
where yt is the number of mobile cellular subscriptions per 100 people at time t, the parameter β1 represents the asymptote toward which the number of subscriptions grows, β2 reflects the relative displacement, and the parameter β3 controls the growth rate (in y 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.Figure2(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.Mobile phone subscriptions (per 100 people) and NLS fit (a), NLS residuals plots (b).
(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β10=100 is therefore reasonable. If we scale time so that t=0 corresponds to 1990, then substituting β10=100 into the model and using the value y0=2.09 from the data, and assuming that the error is 0, we have β20=3.87. Finally, returning to the data, at time t=1 (i.e., in 1991), the number of subscriptions per 100 people was y1=2.95. Using this value along with previously determined start values and again setting the error to 0 gives β30=0.912. The solution to the optimization problem is then reached after six iterations and an achieved convergence tolerance of 4.227×10-6.Results for the naive nonlinear least squares (NLS) estimation of the model are given in Table2. The relative growth rate is estimated at β^3=0.878 and the displacement parameter is estimated at β^2=4.6. More interestingly, the asymptotic upper bound of the S-curve is estimated at β^1=131.76; 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.Table 2
NLS results with errorsut~iid(0,σu2) and MLE results with MA(1) errors given by (35).
Parameters
Naive NLS estimates
ML estimates
Estimates
Stand. errors
Estimates
Stand. errors
β
1
131.8
4.692
129.3
1.543
β
2
4.600
0.135
4.644
0.021
β
3
0.878
0.005
0.8901
0.001
σ
u
2
1.462
0.032
⋯
⋯
γ
⋯
⋯
0.709
0.154
σ
ϵ
2
⋯
⋯
1.074
0.004Figure2(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 ut follows the MA(1) process defined by
(35)ut=ϵt+γϵt-1,withϵt~N(0,σϵ2).Autocorrelation function (a) and partial autocorrelation function (b) of NLS residuals.
(a)
(b)Maximum likelihood (ML) estimates of the model can then be computed by incorporating the above MA(1) 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(1) errors produces different results and smaller standard errors for the estimates. For example, the difference between the NLS estimate for β3 and its ML estimate is 0.0121 and this difference is significant. The moving average coefficient is estimated at γ^=0.709 with a standard error of 0.154. The achieved convergence tolerance is 8.55×10-5.Based on the MLEs, thep-value function p(γ) 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.Table 3
Results for 90% and 95% confidence intervals forγ.
Method
90% CI forγ
95% CI forγ
LR
[0.590, 0.934]
[0.567, 0.978]
Wald
[0.589, 0.887]
[0.546, 0.970]
BN
[0.587, 0.859]
[0.509, 0.911]Figure 4
p-value function for γ.The residuals diagnostics depicted in Figure5 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 -2 and 2, and the ACF plots show that autocorrelation of any lag is not significantly different from zero. This further supports the MA(1) model as a plausible error structure for the data.MLE residuals and their autocorrelation function.
(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
(36)Ct=β1+β2Ytβ3+ut,t=0,1,…,n,
where Yt is personal disposable income and Ct is aggregate personal consumption. The Keynesian consumption function is a restricted version of this model corresponding to β3=1. With this restriction the model is linear. However, if β3 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 β1 and β2 at their linear least squares values that would correspond to the special case of β3=1 and hence use 1 as a starting value for β3. The numerical iterations are therefore started at the linear least squares estimates of β1, β2, and 1 for β3, and convergence is achieved after 8 iterations at a tolerance of 3.111×10-6.Results for the naive linear (OLS) and naive nonlinear least squares (NLS) estimation of the model obtained by assumingut~iid(0,σu2), 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.Table 4
Ordinary least squares, nonlinear least squares, and maximum likelihood with MA(1) errors.
Parameters
OLS estimates
NLS estimates
ML estimates
Estimates
Stand. errors
Estimates
Stand. errors
Estimates
Stand. errors
β
1
−0.3769
0.0914
0.0602
0.6132
0.0775
0.0851
β
2
0.8606
0.0103
0.6986
0.2141
0.6925
0.0082
β
3
1
⋯
1.067
0.0995
1.0701
0.0019
σ
u
2
0.0140
0.0057
0.0145
0.0065
⋯
⋯
γ
⋯
⋯
⋯
⋯
0.1189
0.111
σ
ϵ
2
⋯
⋯
⋯
⋯
0.0123
0.004Except for the intercept the estimated nonlinear relationship gives significant estimates from both the proposed and the traditional approaches (see Table5). In particular, although the power coefficient, β3, 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 0.119, very far from the boundary, implying that the likelihood function is likely well-behaved.Table 5
Results for 95% confidence intervals forβ1, β2, β3, σ2, and γ.
Method
β
1
β
2
β
3
γ
σ
2
LR
[−0.087, 0.234]
[0.672, 0.711]
[1.065, 1.114]
[−0.095, 0.346]
[0.006, 0.023]
Wald
[−0.089, 0.239]
[0.676, 0.708]
[1.066, 1.074]
[−0.098, 0.336]
[0.004, 0.021]
BN
[−0.097, 0.138]
[0.659, 0.721]
[1.056, 1.072]
[−0.019, 0.249]
[0.010, 0.016]However, the residuals do not appear to display serial correlation as suggested by Figure6. 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(1) structure, as evident by the confidence intervals obtained in Table 5.NLS residuals and their autocorrelation function.
(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.


## 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
(32)yt=θ1[1+eθ2+θ3xt]-1+ut,t=1,…,n
with disturbances ut generated by the MA(1) process defined by
(33)ut=ϵt+γϵt-1,t=…,-1,0,1,….The errorsϵt are assumed to be distributed independently as standard normal. The explanatory variable xt is a vector of deterministic integers running from 0 to n-1. The parameter values are given as θ1=56, θ2=2.9, and θ3=-0.24. Samples of sizes n=10 and n=20 are used along with 5,000 replications each. The moving-average coefficient takes values γ∈{-0.9;-0.6;-0.3;0;0.3;0.6;0.9} and confidence intervals of 90% and 95% are computed for each value of the moving average parameter. Note that the value γ=0 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. Table1 reports simulation results for a sample size of n=10, while Figure 1 graphically illustrates the results when the sample size is n=20.Table 1
Simulation results forn=10, 95% CI.
γ
Method
Central  coverage
Upper  error
Lower  error
Average  bias
−0.9
LR
0.9120
0.0685
0.0195
0.0245
Wald
0.9155
0.0405
0.0440
0.0173
BN
0.9385
0.0285
0.0330
0.0058
−0.6
LR
0.9285
0.0655
0.0060
0.0298
Wald
0.9250
0.0560
0.0190
0.0185
BN
0.9390
0.0275
0.0365
0.0070
−0.3
LR
0.9180
0.0610
0.0210
0.0200
Wald
0.9015
0.0690
0.0295
0.0243
BN
0.9405
0.0260
0.0335
0.0047
0
LR
0.8510
0.0585
0.0905
0.0495
Wald
0.8150
0.0805
0.1045
0.0675
BN
0.9365
0.0280
0.0355
0.0068
0.3
LR
0.8650
0.0760
0.0590
0.0425
Wald
0.8230
0.0825
0.0945
0.0635
BN
0.9380
0.0295
0.0325
0.0060
0.6
LR
0.9090
0.0370
0.0540
0.0205
Wald
0.8730
0.0410
0.0860
0.0385
BN
0.9360
0.0305
0.0335
0.0070
0.9
LR
0.8890
0.0170
0.0940
0.0385
Wald
0.8630
0.0215
0.1155
0.0470
BN
0.9345
0.0295
0.0370
0.0083Figure 1
Graphical illustration of the simulation results forn=20, 95% CI.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 Figure1 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
(34)yt=β1exp⁡[-β2β3t]+ut,t=0,1,…,n,
where yt is the number of mobile cellular subscriptions per 100 people at time t, the parameter β1 represents the asymptote toward which the number of subscriptions grows, β2 reflects the relative displacement, and the parameter β3 controls the growth rate (in y 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.Figure2(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.Mobile phone subscriptions (per 100 people) and NLS fit (a), NLS residuals plots (b).
(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β10=100 is therefore reasonable. If we scale time so that t=0 corresponds to 1990, then substituting β10=100 into the model and using the value y0=2.09 from the data, and assuming that the error is 0, we have β20=3.87. Finally, returning to the data, at time t=1 (i.e., in 1991), the number of subscriptions per 100 people was y1=2.95. Using this value along with previously determined start values and again setting the error to 0 gives β30=0.912. The solution to the optimization problem is then reached after six iterations and an achieved convergence tolerance of 4.227×10-6.Results for the naive nonlinear least squares (NLS) estimation of the model are given in Table2. The relative growth rate is estimated at β^3=0.878 and the displacement parameter is estimated at β^2=4.6. More interestingly, the asymptotic upper bound of the S-curve is estimated at β^1=131.76; 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.Table 2
NLS results with errorsut~iid(0,σu2) and MLE results with MA(1) errors given by (35).
Parameters
Naive NLS estimates
ML estimates
Estimates
Stand. errors
Estimates
Stand. errors
β
1
131.8
4.692
129.3
1.543
β
2
4.600
0.135
4.644
0.021
β
3
0.878
0.005
0.8901
0.001
σ
u
2
1.462
0.032
⋯
⋯
γ
⋯
⋯
0.709
0.154
σ
ϵ
2
⋯
⋯
1.074
0.004Figure2(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 ut follows the MA(1) process defined by
(35)ut=ϵt+γϵt-1,withϵt~N(0,σϵ2).Autocorrelation function (a) and partial autocorrelation function (b) of NLS residuals.
(a)
(b)Maximum likelihood (ML) estimates of the model can then be computed by incorporating the above MA(1) 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(1) errors produces different results and smaller standard errors for the estimates. For example, the difference between the NLS estimate for β3 and its ML estimate is 0.0121 and this difference is significant. The moving average coefficient is estimated at γ^=0.709 with a standard error of 0.154. The achieved convergence tolerance is 8.55×10-5.Based on the MLEs, thep-value function p(γ) 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.Table 3
Results for 90% and 95% confidence intervals forγ.
Method
90% CI forγ
95% CI forγ
LR
[0.590, 0.934]
[0.567, 0.978]
Wald
[0.589, 0.887]
[0.546, 0.970]
BN
[0.587, 0.859]
[0.509, 0.911]Figure 4
p-value function for γ.The residuals diagnostics depicted in Figure5 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 -2 and 2, and the ACF plots show that autocorrelation of any lag is not significantly different from zero. This further supports the MA(1) model as a plausible error structure for the data.MLE residuals and their autocorrelation function.
(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
(36)Ct=β1+β2Ytβ3+ut,t=0,1,…,n,
where Yt is personal disposable income and Ct is aggregate personal consumption. The Keynesian consumption function is a restricted version of this model corresponding to β3=1. With this restriction the model is linear. However, if β3 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 β1 and β2 at their linear least squares values that would correspond to the special case of β3=1 and hence use 1 as a starting value for β3. The numerical iterations are therefore started at the linear least squares estimates of β1, β2, and 1 for β3, and convergence is achieved after 8 iterations at a tolerance of 3.111×10-6.Results for the naive linear (OLS) and naive nonlinear least squares (NLS) estimation of the model obtained by assumingut~iid(0,σu2), 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.Table 4
Ordinary least squares, nonlinear least squares, and maximum likelihood with MA(1) errors.
Parameters
OLS estimates
NLS estimates
ML estimates
Estimates
Stand. errors
Estimates
Stand. errors
Estimates
Stand. errors
β
1
−0.3769
0.0914
0.0602
0.6132
0.0775
0.0851
β
2
0.8606
0.0103
0.6986
0.2141
0.6925
0.0082
β
3
1
⋯
1.067
0.0995
1.0701
0.0019
σ
u
2
0.0140
0.0057
0.0145
0.0065
⋯
⋯
γ
⋯
⋯
⋯
⋯
0.1189
0.111
σ
ϵ
2
⋯
⋯
⋯
⋯
0.0123
0.004Except for the intercept the estimated nonlinear relationship gives significant estimates from both the proposed and the traditional approaches (see Table5). In particular, although the power coefficient, β3, 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 0.119, very far from the boundary, implying that the likelihood function is likely well-behaved.Table 5
Results for 95% confidence intervals forβ1, β2, β3, σ2, and γ.
Method
β
1
β
2
β
3
γ
σ
2
LR
[−0.087, 0.234]
[0.672, 0.711]
[1.065, 1.114]
[−0.095, 0.346]
[0.006, 0.023]
Wald
[−0.089, 0.239]
[0.676, 0.708]
[1.066, 1.074]
[−0.098, 0.336]
[0.004, 0.021]
BN
[−0.097, 0.138]
[0.659, 0.721]
[1.056, 1.072]
[−0.019, 0.249]
[0.010, 0.016]However, the residuals do not appear to display serial correlation as suggested by Figure6. 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(1) structure, as evident by the confidence intervals obtained in Table 5.NLS residuals and their autocorrelation function.
(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 accuratep-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(1) 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(1) 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.



---

*Source: 207087-2014-02-13.xml*