Two-Stage Method Based on Local Polynomial Fitting for a Linear Heteroscedastic Regression Model and Its Application in Economics

We introduce the extension of local polynomial fitting to the linear heteroscedastic regression model. Firstly, the local polynomial fitting is applied to estimate heteroscedastic function, then the coefficients of regression model are obtained by using generalized least squares method. One noteworthy feature of our approach is that we avoid the testing for heteroscedasticity by improving the traditional two-stage method. Due to nonparametric technique of local polynomial estimation, we do not need to know the heteroscedastic function. Therefore, we can improve the estimation precision, when the heteroscedastic function is unknown. Furthermore, we focus on comparison of parameters and reach an optimal fitting. Besides, we verify the asymptotic normality of parameters based on numerical simulations. Finally, this approach is applied to a case of economics, and it indicates that our method is surely effective in finite-sample situations.


Introduction
The heteroscedasticity in classical linear regression model is defined by the variances of random items which are not the same for different explanatory variables and observations.Heteroscedasticity often occurs in data sets in which there is a wide disparity between the largest and smallest observed values.The larger the disparity between the size of observations in a sample, the larger the likelihood that the error term observations associated with them will have different variances and therefore be heteroscedastic.That is, we would expect that the error term distribution for very large observations might have a large variance, but the error term distribution for small observations might have a small variance.Besides, researchers have observed that heteroscedasticity is usually found in cross-sectional data rather than in time series data 1, 2 .In cross-sectional data we generally deal with members of a population at a given point in time, such as individual consumers or their families; firms; industries; or geographical subdivisions, such as small, medium, or large firms, or low, medium, or high income.Besides, this study is also more common in economics, such as the relationship of sale and research and development, and profit in a certain year.
When there is a heteroscedasticity in a linear regression model 3-6 , we can apply parametric methods, also the heteroscedastic approaches can be applied, such as Park test, White test and so on.Another heteroscedastic method is shown in 7 .When parametric approaches are applied, estimations of parameters we obtained by ordinary least squares estimation OLS are still linear and unbiased.However, the efficiency is bad 8-10 .This could lead to a uncorrect statistical diagnosis for the parameters' significance test.Similarly, it unecessarily enlarges the confidence interval when we estimate the parameter interval.Besides, the accuracy of predictive value may lower when we estimate with the regression model that we obtained.In order to solve the problem above, we can use generalized least squares estimation GLS when the covariance matrix of the random items is known.If it is unknown, we usually use two-stage least squares estimate, that is, we first estimate variances of the residual error, and then the generalized least squares estimator is used to obtain the coefficients of the model by using the estimate of variances of the random items 11, 12 .However, the traditional estimation method is that we suppose the residual error variances as a certain parametric model.In this paper, we try to applying local polynomial fitting to random item variances as the first step, and then GLS is used to estimate the coefficients of the model.A problem cannot be neglected is that the existence of inverse of regressor matrix transpose multiplied by regressor matrix is necessary for the existence of the unique estimator when applying the least squares method or their generalized versions.In the generalized versions, because kernel function is a symmetric probability density function with bounded support, the weighted matrix based on kernel function is positive definite, generally being also symmetric as done in usual theoretical and practical selections.In other words, the regressor matrix has to be full column rank and weighted matrix has to be nonsingular.The issues have been discussed in many papers used in theoretical issues and many applications, see 13-17 .On the one hand, because of local polynomial fitting's various nice statistical properties, the estimations obtained with this technology also possess the same good statistical properties 12, 18, 19 .On the other hand, we exploit a heteroscedastic regression model rather than the artificial structure of heteroscedasticity.Then, we can directly get the heteroscedastic function based on the nonparametric technique, which shows the relationship between variance function of random items and explanatory variables from regression results.Thus, it is unnecessary to test heteroscedasticity of the model.Particularly, the estimated value by local polynomial fitting is more accurate than that by the traditional method.Besides, we study variance function fitting when parameters change and reach the optimal parameters.Finally, a case of economics is cited in order to show that our method is indeed effective in finitesample situations.
The rest of this paper is organized as follows.In Section 2 we construct local polynomial regression: in Section 2.1 we talk about local polynomial fitting, in Section 2.2 we study estimation with Parameters selections.Section 3 contains two-stage method with local polynomial fitting.In Section 4, we do the simulations on a given model and study the fittings under different parameters.In Section 5, we collect some real data and use the local polynomial estimating the coefficients.We draw the results in Section 6.

Local Polynomial Regression
Local polynomial regression 20-23 is a widely used nonparametric technique.Local polynomial fitting is an attractive method both from theoretical and practical point of view.
Local polynomial method has a small mean squared error compared with the Nadaraya-Watson estimator which leads to an undesirable form of the bias and the Gasser-Muller estimator which has to pay a price in variance when dealing with a random design model.Local polynomial fitting also has other advantages.The method adapts to various types of designs such as random and fixed designs, highly clustered and nearly uniform designs.Furthermore, there is an absence of boundary effects: the bias at the boundary stays automatically of the same order as the interior, without use of specific boundary kernels.The local polynomial approximation approach is appealing on general scientific grounds: the least squares principle to be applied opens the way to a wealth of statistical knowledge and thus easy generalizations.In this section, we briefly outline the idea of the extension of local polynomial fitting to linear regression.

Local Polynomial Fitting
Consider the bivariate data X 1 , Y 1 , . . ., X n , Y n , which form an independent and identically distributed sample from a population X, Y .Of interest is to estimate the regression function m x 0 E Y | X x 0 and its derivatives m x 0 , m x 0 , . . ., m p x 0 .To help us understand the estimation methodology, we can regard the data as being generated from the model where E ε 0, Var ε 1, and X and ε are independent.However, this location-scale model assumption is not necessary for our development, but is adopted to provide intuitive understanding.We always denote the conditional variance of Y given X x 0 by σ 2 x 0 and the marginal density of X, that is, the design density, by f • .
Suppose that the p 1 th derivative of m x at the point x 0 exists.We then approximate the unknown regression function m x locally by a polynomial of order p.A Taylor expansion gives, for x in a neighborhood of x 0 , This polynomial is fitted locally by a weighted least squares regression problem: minimize where h is a bandwidth controlling the size of the local neighborhood, and K h • K •/h /h with K a kernel function assigning weights to each datum point.Throughout this paper we assume that K is a symmetric probability density function with bounded support, although this technical assumption can be relaxed significantly 24-26 .Denote by β j , j 0, . . ., p, the solution to the least squares problem 2.3 .It is clear from the Taylor expansion in 2.2 that m ν x 0 ν! β ν is an estimator for m ν x 0 , ν 0, 1, . . ., p.To estimate the entire function m ν • we solve the above weighted least squares problem for all points x 0 in the domain of interest.
It is more convenient to work with matrix notation.Denote by X the design matrix of problem 2.3 : and let Further, let W be the n × n diagonal matrix of weights: Then the weighted least squares problem 2.3 can be written as min with β β 0 , . . ., β p T .The solution vector is provided by weighted least squares theory and is given by where the regressor matrix X has to be full column rank and the weighted matrix W has to be nonsingular.If such an inverse in 2.8 does not exist, a method of ridge regression can be adopted to solve this problem 27 .Furthermore, we can get the estimation m x 0 , where E 1 is a column vector the same size of β with the first element equal to 1, and the rest equal to zero, that is, E 1 1, 0, . . ., 0 1× p 1 .
Computing the β will suffer from large computational cost.We can use the recursive least squared method to reduce the computation complexity, and it is very powerful especially in the local polynomial fitting problems.There are several important issues about the bandwidth, the order of local polynomial function and the kernel function which have to be discussed.The three problems will be presented in Section 2.2.

Parameters Selections
To implement the local polynomial estimator, one need to choose the order p, the kernel K and the bandwidth h.These parameters are of course related to each other.
First of all, the choice of the bandwidth parameter h is considered, which plays a rather crucial role.A too large bandwidth under-parametrizes the regression function 28, 29 , causing a large modeling bias, while a too small bandwidth over-parametrizes the unknown function and results in noisy estimates.The basic idea is to find a bandwidth h that minimizes the estimated mean integrated square error MISE : where X X 1 , . . ., X p , the asymptotic bias and the asymptotic variance are denoted in Lemma A.1 of the appendix.Then we can find an asymptotically optimal constant bandwidth given by where, c ν,p K is a constant which relates to the kernel function and the order of the local polynomial, f x , is the density function of m • .However, this ideal bandwidth is not directly usable since it depends on unknown functions.Another issue in local polynomial fitting is the choice of the order of the polynomial.Since the modeling bias is primarily controlled by the bandwidth, the issue is less crucial however.For a given bandwidth, h, a large value of p would expectedly reduce the modeling bias, but would cause a large variance and a considerable computational cost.It is shown in 30 that there is a general pattern of increasing variability: for estimate m ν x 0 , there is no increase in variability when passing from an even i.e., p − ν even p ν 2q order fit to an odd p ν 2q 1 order fit, but when passing from an odd p ν 2q 1 order fit to the consecutive even p ν 2q 2 order fit, there is a price to be paid in terms of increased variability.Therefore, even order fits p ν 2q are not recommended.Since the bandwidth is used to control the modeling complexity, we recommend the use of lowest odd order, that is, p ν 1, or occasionally p ν 3.
Another question concerns the choice of the kernel function K. Since the estimation is based on the local regression 2.3 , no negative weight K should be used.As shown in 30 , the optimal weight function is K z 3/4 1 − z 2 , the Epanechnikov kernel, which minimizes the asymptotic mean square error MSE of the resulting local polynomial estimators.

Two-Stage Method with Local Polynomial Fitting
Let the dependent variable x and the explanatory variable y fulfill the following regression model: where y i are the observations and x i , i 1, . . ., n are independent variables.Denote Therefore, 3.1 can be abbreviated as y X x β σ x ε.

3.3
Suppose that are not all equal, that is, there is a heteroscedasticity in model 3.4 .Therefore, GLS for β is If covariance matrix Σ is known, the coefficients β can be estimated.Equation 3.4 is considered as the weighted least squares estimation WLS for β and it possesses nice qualities.However, how to estimate the Σ is still a problem.Therefore, the so-called two-stage method of estimation is used to solve the heteroscedasticity problem.The two-stage method based on local polynomial fitting can be depicted as follows: firstly, apply local polynomial fitting to get the estimate for , that is, σ 2 i for σ 2 i , and then we can obtain the estimate β for β by using 3.5 .The estimator follows that we construct the following regression model in order to estimate σ 2 i , where u i is the difference between σ x ε i 2 and its expectation.Suppose that b is the OLS of model 3.4 .Although the ordinary least squares estimate b is ineffective, it is still consistent.Therefore, the corresponding residuals hold that

3.7
Consequently, we can approximately get

3.8
It can be taken as a regression model, in which the variance function is regression function and the squared residuals e 2 i are dependent variables.In order to estimate this model, parameter estimation method would usually be taken in some articles.In other words, they suppose σ 2 i f c T , x i , where the form of f is known and c c 0 , c 1 T are the parameters to be estimated.Note that what we usually discuss about are σ 2 σ 2 c T , x i , σ i 2 σ 2 i ln c T , x i and so on 31 .However, the discussion for these models requires the analysts to have a better understanding of the background in practical problems.As an example, variance of asset return is always in direct proportion with its past return.Since the variance function must be nonnegative, a nonparametric method is proposed to fit σ 2 i .This method can be depicted as follows.Then, a p-order local polynomial estimation for the variances function σ 2 x is obtained according to formula 2.2 .Using the least squares method for the data around the local window, we can estimate the local intercept via minimizing with respect to {α i } p i 0 .Therefore, the solution vector can be written as

3.10
where the design matrix the weighted matrix 12 and e 2 e 2 1 , e 2 2 , . . ., e 2 n T .Consequently, the estimated variance function is σ 2 x E 1 α.Finally, we can get two-stage estimate β for β by substituting estimate σ 2 i for σ 2 i into 3.5 , which has some wonderful statistics qualities, seeing Lemma A.2, Lemma A.3, and Theorem A.4 in the appendix.

Simulation and Analysis
In order to discuss the qualities of β under the limited sample, this section gives the following model, in which we do the comparison and study the fittings under different parameters.Considering the practical background which is applied to economics, we suppose the variance function of the following form.Besides, we study the fitting effects of variance in different parameters to obtain the best one.
Denote the linear model by where x i i 1, . . ., n are independent variables, y i i 1, . . ., n are observations.Also, E σ x i ε i σ x i ε i T diag σ 2 1 , . . ., σ 2 n , and σ 2 1 , . . ., σ 2 n are not all equal.Suppose that the variance function of the error term is σ 2 x e 0.6x .
Step 1. Firstly, obtain the estimation of two coefficients in model 4.1 with the ordinary least squares estimation.Secondly, calculate the squares of the residuals e 2 e 2 1 , . . ., e 2 n T .Thirdly, do the local polynomial regression based on the model 3.8 .In this section, it is necessary to discuss how to choose parameters.Suppose that the range of x is −1, 5 and the kernel function is the Epanechnikov kernel K z 3/4 1−z 2 .In addition, the criteria of selecting  a bandwidth h is minimizing the mean integrated square error MISE , which can be given by We could get values of MISE under different orders and bandwidths by 10000 replicates calculation with the above equation, where n 300, see Table 1.
Furthermore, the scatter plot about h and MISE under different orders can be drawn, see Figure 1.
From Figure 1, it can be seen that values of MISE are the minimum when h 0.5 whether p 1, p 2, and p 3. Therefore, h 0.5 is the optimal bandwidth.Further, value of MISE when p 2 is the minimum among the above three, that is, p 2 is the optimal order.It can be drawn the conclusion that the optimal parameters are h 0.5 and p 2, respectively.Figure 2 shows the fitting plot and residual plot for variance function after 10000 replicates, where h 0.5, p 2 and n 300.
Step 2. Now we substitute σ 2 i which is obtained from Step 1 into model 3.5 , and then we will get the estimation β of GLS, that is, β 0 3.1004, β 1 1.4987.Figure 3 depicts the histograms and asymptotic distributions of β 0 and β 1 , by which we do 10000 replicates with GLS and  choose n 300.It is easy to see from Figure 3 that the estimated distributions of parameters are subject to normal asymptotically.Besides, the OLS for β 0 and β 1 can easily be obtained, say, b 0 3.1813, b 1 1.4394.The fitted and true values for GLS and ones of OLS are listed in Table 1.Here, the relative error is defined by Re | ξ − ξ|/ξ × 100%, where ξ and ξ present true and fitted values, respectively.From the comparison in Table 2, we can conclude that the estimators for parameters by GLS are much better than those by OLS.Furthermore, curves of original and regression functions estimated by GLS and OLS are plotted together in order to demonstrate accuracy by GLS, see Figure 4.It is not difficult to see that the GLS regression curve is almost coincident with original curve, while there is a parent bias between OLS regression and original curve.Consequently, estimators for parameters by GLS are what we require.

Application
As an example of pure cross-sectional data with potential for heteroscedasticity, consider the data given in Table 3, which gives data on per capita consumption Y and per capita gross domestic product GDP X for 31 province or city in the People's Republic of China in 2008.Since the cross-sectional data presented in this table are quite heterogenous in a regression of per capita consumption Y on per capita GDP X , heteroscedasticity is likely.
If we want to understand relationship between per capita consumption and per capita GDP, then the regression function is as follows: where B 0 and B 1 are regression coefficients and μ is a vector of random errors.The ordinary least squares OLS estimations for B 0 and B 1 can easily be obtained, that is, B 0 679.72 and B 1 0.31117.As follow, according to the known sample data, the regression plot with OLS can easily be drawn, seeing Figure 5.Not surprising, there is a positive relationship between per capita consumption and per capita GDP, although it is not statistically significant at the traditional levels.To testify if the regression 5.1 suffers from heteroscedasticity, we obtain the residuals of the model and plotted them against with per capita GDP, as shown in Figure 6, in which we can see that the residuals are increasing around horizontal axis with the increase in per capita GDP.Therefore, there is a heteroscedasticity in the regression 5.1 .Heteroscedasticity can be also obtained through other tests such as Park test and White test.So it can be said that the regression 5.1 suffers from heteroscedasticity.Furthermore, the generalized least squares GLS estimations B 0 665.84 and B 1 0.31265 can easily be got after 10000 replicates,   provided with the bandwidth h 0.5, the order p 2, and MISE 0.275746.Then the regression plot can be drawn in Figure 7. Compared with Figures 5 and 7, it can be said that the distribution of points around regression line in Figure 7 is more uniform than that in Figure 5. Finally, it can be exactly said that more accurate regression can be obtained by GLS than that by OLS.

Conclusions
In this paper we presented a new method for estimation of linear heteroscedastic regression model based on local polynomial estimation with nonparametric technique.The proposed scheme firstly adopted the local polynomial fitting to estimate heteroscedastic function, then the coefficients of regression model are obtained based on generalized least squares method.Our approach avoided the test of heteroscedasticity for the linear model.Due to nonparametric technique of local polynomial estimation, if the heteroscedastic function is unknown, the precision of estimation was improved.Furthermore, the effect of parameters on the fitting was researched and the optimal fitting was obtained.Besides, the asymptotic normality of parameters was verified by the results of numerical simulations.Finally, the simulation results under different parameters and local polynomial estimation of real data in a case of economics really indicated that our approach was effective in finite-sample situations, which did not need to assume the form of heteroscedastic function.The presented algorithm could be easily used to heteroscedastic regression model in some practical problems.Lemma A.2.For the model 3.4 , ε is assumed to be a vector and follows normal distribution with mean zero and covariance Σ, where Σ > 0, and P lim n → ∞ X Σ −1 X/n Q, a finite positive definite matrix.Then one has 1 β is a consistent estimator of β.

The proof of Lemma
2 β is asymptotically normal with the mean β and covariance matrix X Σ −1 X −1 .
The proof is shown in 32 .
Lemma A.3.Suppose σ 2 x is a p-order local polynomial estimator (LPE) of the variance function σ 2 x .Under conditions h n → 0, and nh n → ∞: 1/nh n .
The proof is shown in 18 .
Theorem A.4.For the model 3.4 , ε is assumed to be a vector and follows normal distribution with mean zero and covariance Σ, where Σ > 0, and p lim n → ∞ X T X/n Q 1 , p lim n → ∞ X T ε/ √ n Q 2 , where Q 1 and Q 2 are positive definite matrixes.Besides, they satisfy the conditions of Lemma A.3, then β LPE is asymptotically normal with mean β and covariance X T Σ −1 X −1 .
Proof.By Lemma A.2 and the bibliography 32 , in order to prove the above results, we only need to testify where Σ is an estimator of Σ, that is to testify A.9 By Lemma A.3, whether p is an odd or an even, we have The conclusions follow.

Figure 2 :
Figure 2: The fitting plot and its residual plot when h 0.5 and p 2: a fitting plot, b residual plot.

Figure 4 :
Figure 4: Curves of original and regression functions estimated by GLS and OLS.
A.1 is shown in 23 .An immediate consequence of Lemma A.1 is that the asymptotic bias and the asymptotic variance for the local polynomial estimator are defined as mean square error MSE at point x is MSE x AB x 2 AV x .A.4

Table 1 :
Values of MISE under different orders and bandwidths.

Table 2 :
Two-stage estimates for β 0 and β 1 of GLS and OLS.

Table 3 :
Per capita consumption and per capita GDP for all provinces or cities in the People's Republic of China, 2008 values are in yuan of RMB .
an odd, and the variance function σ 2 x is continuous p 1 derivative and bounded, then if p is an even, and the variance function σ 2 x is continuous p 1 and p 2 derivative and bounded, then where K 1 , K 2 , K 3 , and K 4 are constants.Then we can get max 1≤i≤n MSE σ 2