Instrumental Variable Quantile Regression of Spatial Dynamic Durbin Panel Data Model with Fixed Effects

: This paper studies a quantile regression spatial dynamic Durbin panel data (SDDPD) model with ﬁxed effects. Conventional ﬁxed effects estimators of quantile regression speciﬁcation are usually biased in the presentation of lagged response variables in spatial and time as regressors. To reduce this bias, we propose the instrumental variable quantile regression (IVQR) estimator with lagged covariates in spatial and time as instruments. Under some regular conditions, the consistency and asymptotic normalityof the estimators are derived. Monte Carlo simulations show that our estimators not only perform well in ﬁnite sample cases at different quantiles but also have robustness for different spatial weights matrices and for different disturbance term distributions. The proposed method is used to analyze the inﬂuencing factors of international tourism foreign exchange earnings of 31 provinces in China from 2011 to 2017.


Introduction
A panel data set is one that follows a given sample of individuals over time and thus provides multiple observations on each individual in the sample (Hsiao, 2014 [1]). Compared with cross-sectional or time series data models, panel data models contain more information and variation and have greater capacity for dealing with complex situations. The specification and estimation of panel data models have received substantial attention since the seminal paper (Balestra and Nerlove [2]) was published. Their theories, methods, and applications of panel data models are abundant and can be found in books by Baltagi [3], Arellano [4], Baldev and Baltagi [5], and Hsiao [1], among others. Adding a lagged response variableinto a panel data model, we can naturally establish a dynamic panel data model that is more flexible. The estimation of dynamic panel data models has been considered in many studies. Nickell [6] pointed out that ordinary least squares estimators in dynamic panel data models are biased and inconsistent. To overcome the adverse effect of lagged response variable, statisticians have developed the following three commonly used estimation approaches to study the dynamic panel data models: maximum likelihood (ML) estimation (Nerlove [7]; Hsiao et al. [8]), generalized method of moment (GMM) estimation (Arellano and Bond [9]; Blundell and Bond [10]), and least square dummy variable (LSDV) estimation (Kiviet [11]; Hahn and Kuersteiner [12]).
By using spatial weights matrix (for details, see Anselin [13] and Elhorst [14]) to describe the presence and strength of a link between observations in the sample, spatial econometrics deals with the spatial correlations between different individual units in panel data models. The spatial correlations of observations have been widely used to describe economic activities such as spillover effects, social interactions, externalities, etc. According to Elhorst [14], there are three basic categories of spatial correlations that have been commonly used in the literature. The first category is endogenous interaction effects among observations of response variable at different individual units, where the observations of response variable for one individual unit is jointly determined with that of neighboring units.By allowing endogenous interaction effects in a dynamic panel data model, Anselin [13] provided a framework for spatial dynamic panel data (SDPD) models. Yu et al. [15] studied the asymptotic properties of quasi-maximum likelihood (QML) estimators for a SDPD model with individual fixed effects. Lee and Yu [16] extended this methodology to a SDPD model with individual and time fixed effects. Yu and Lee [17] obtained the asymptotic properties of QML estimators for a unit root SDPD model. Lee and Yu [18] studied the asymptotic properties of GMM estimators for a SDPD model. The second category is exogenous interaction effects among observations of covariates at different individual units, where the response variable of one individual unit depends on independent covariates of other units.By considering spatially lagged contemporaneous covariates and spatially time-lagged covariates, which are the so-called spatial Durbin terms (see LeSage and Pace [19]), we can capture the direct externality effects of exogenous variables on response variable. There are many literature on the applications of dynamic spatial Durbin model (DSDM); see Basile et al. [20], Vega and Elhorst [21], Rios [22], and Feng and Wang [23], among others. The third category is interaction effects among disturbance terms of different individual units. That is, the determinants of the response variable omitted from the model are spatially autocorrelated or with a situation where unobserved shocks follow a spatial pattern (Elhorst [14], 2014). The allowance of random error spatial dependence complicates dynamic panel data models. Elhorst [24], Mutl [25], Yang et al. [26] and Su and Yang [27] studied the dynamic panel data models with spatial errors. By a suitable combination of the above three basic categories of spatial interaction effects, some authors established different types of SDPD models and obtained their estimation methods; see Lee and Yu [28], Qu et al. [29], Shi and Lee [30], and Yang and Lee [31], among others.
Quantile regression was first proposed by Koenker and Bassett [32]; it is an extension of mean regression to the quantiles of responses. It offers a systematic strategy for describing how covariates affect the location, scale, and shape of response distribution. Compared with mean regression, which explores how the average of response variable dependson the values of conditioning covariates, quantile regression has several advantages. First, quantile regression allows for separate modelling at different quantiles of response distribution such that the impacts of covariates can be differentiated at different quantiles of response distribution. As a result, we can have a much richer view in applications than could be achieved by mean regression. Second, quantile regression does not need any reliance on global distributional assumptions, which is an attractive feature of quantile regression compared with mean regression. Third, quantile regression can fully describe the heterogeneity in response as well as takes unobserved heterogeneity into consideration. Forth, quantile estimators are robust and less sensitive to the heteroscedasticity and outliers that frequently present in real data. Based on the above advantages, quantile regression has gradually become a valuable statistical methodology and attracted considerable interest in many research fields. Koenker [33] provided an excellent exposition of quantile regression.
Quantile regression for panel data models can not only fully describe the conditional distribution of a response variable but also control individual specific heterogeneity via fixed effects. Koenker [34] explored quantile regression for longitudinal data with individual fixed effects. The fixed effects were regarded as location shift parameters of conditional quantiles and shrank toward a common value. Based on 1 regularization, the author proposed a penalized quantile regression estimator and obtained its asymptotic properties. Galvao [35] extended quantile regression to dynamic panel data model. IVQR, which was proposed by Chernozhukov and Hansen [36], was selected to reduce the estimation bias caused by dynamic term. Based on the estimated model, the author obtained its consistent estimator, asymptotic properties, and prediction. The estimation technique was used to predict the output growth rates of 18 OECD countries. Tao et al. [37] considered Hausman-Taylor instrumental variables and obtained two quantile regression estimators for a dynamic panel data model, the asymptotic properties of the proposed estimators were studied, and the influencing factors of commercial residential prices in large and moderate cities in China were analyzed. To the best of our knowledge, only two papers on quantile regression of the spatial panel data models have been published. Dai et al. [38] studied the IVQR of a spatial error panel data model with individual fixed effects. Dai et al. [39] extended the IVQR to a spatial autoregressive panel data model with autoregressive disturbances. The asymptotic properties of estimators in these two papers are based on the assumption that the observations of response variable are independent and identically distributed. However, in the spatial panel data models, the assumptions are whether the endogenous interaction in spatial dimensions or autocorrelation in time series makes the independent among observations of response variable unreasonable.
To summarize, our paper provides a complete set of quantile regression estimation methodologies for the SDDPD model with fixed effects, accommodating different correlations among response variable and covariates (spatial, temporally, and spatiotemporal). We find that the panel data fixed effects estimators of quantile regression specification are usually biased. To reduce the bias, we suggest the use of the IVQR method by Chernozhukov and Hansen [36] along with lagged covariates in spatial and time as instruments. Thus, the IVQR estimator combines the instrument variable concept for a SDDPD model and quantile regression framework. We study the asymptotic properties of quantile regression estimators under reasonable assumptions. Monte Carlo simulations under various quantiles indicate that our estimators perform well in finite samples even if the random error distribution is heavy-tailed or asymmetrical. Finally, we illustrate our method with an application to analyze the factors affecting international tourism foreign exchange earnings of 31 provinces in China from 2011 to 2017. Our method is quite general and can be generalized to other types of SDPD models and DSDM models. It is relatively easy to be used in analyzing real data sets and thus greatly facilitates empirical research. This paper is organized as follows. Section 2 introduces the IVQR of SDDPD models with fixed effects, and its estimators are constructed. Section 3 proves the asymptotic properties of the IVQR estimators under some regular conditions. Section 4 reports the finite sample performance of estimates by Monte Carlo simulations. Section 5 illustrates the proposed method by exploring the influencing factors of international tourism foreign exchange earnings in 31 provinces of China from 2011 to 2017.

Model and Estimation
Consider a spatial dynamic Durbin panel data (SDDPD) model with fixed effects under quantile restrictions: where y t = (y 1t , · · · , y Nt ) is an N × 1 vector consisting of the observed values of the cross sectional individuals (i = 1, · · · , N) at time t (t = 1, · · · , T), X t = (x 1t , · · · , x Nt ) is an N × p matrix of covariates at time t, η 0 = (η 01 , · · · , η 0N ) is an N × 1 vector consisting of individual fixed effects, and ε t = (ε 1t , · · · , ε Nt ) is an N × 1 vector of innovation values at time t. Normalized spatial weights matrix W N is an N × N known constant matrix, its diagonal elements are all 0, and W N times a vector or matrix denotes its spatially lagged value. Serially lagged value is represented by the subscript t − 1 of the vector or matrix. The scale parameters λ(τ)(|λ(τ)| < 1), γ 1 (τ)(|γ 1 (τ)| < 1), and γ 2 (τ)(|γ 2 (τ)| < 1) are the parameters of the response variable lagged in space, W N y t ; the response variable lagged in time, y t−1 ; and the response variable lagged in space and time, W N y t−1 , respectively. The p × 1 vectors β 1 (τ), β 2 (τ), β 3 (τ) and β 4 (τ) are parameters of exogenous variables. All parameters depend on τ. This model can be used to simultaneously study the spatial correlation and dynamic effects of response variable and covariates, and the spatial correlation of the lagged response variable and covariates. Elhorst [14] pointed out its application values.
Extending the set of response variable to include endogenous interaction effects, y t−1 (y t−1 , W N y t−1 ), and the set of covariates to include exogenous interaction effects, X t (X t , W N X t , X t−1 , W N X t−1 ). Let γ(τ) γ 1 (τ), γ 2 (τ) and β(τ) β 1 (τ), β 2 (τ), β 3 (τ), β 4 (τ) ; then, the reduced form of model (1) is as follows: The model (2) can be written in detail as where y it is an observation of the response variable for subject i at time t, x it is a p × 1 vector composed of observations of exogenous variables, η 0i is the fixed effect of the ith individual, and ε it is a random error term. The lagged values of the response variable y it in space and time are represented by ∑ N j=1 w ij y jt and y i,t−1 respectively, where w ij denotes the (i, j)th element of spatial weights matrix W N . The scale parameters λ(τ) (|λ(τ)| < 1) and γ(τ) (|γ(τ)| < 1) characterize the spatial effect and dynamic effect, respectively, and β(τ) is a p × 1 regression coefficient vector. We always denote that y it = 0 when t ≤ 0.
Let y = (y 1 , · · · , y T ) be an NT × 1 vector, X = (X 1 , · · · , X T ) be an NT × p matrix, and ε = (ε 1 , · · · , ε T ) be an NT × 1 vector; then, the model (3) can represented as a more concise matrix form where W = I T W N is an NT × NT weights matrix, L is the lagged operator, M = l T I N is an NT × N matrix, I k denotes a k × k identity matrix, l T denotes a T × 1 vector with each element is equal to 1, and represents the Kronecker product. Denote H(λ, γ) = I NT − λW − γL; then, (4) can be reduced as follows: provided that H λ(τ), γ(τ) is nonsingular. This form is used to derive the asymptotic properties of estimators proposed below. It is obvious that the models (2)-(5) are equivalent, and the model (1) can be transformed into the model (2). Therefore, to study the quantile regression model (1), we just need to study the model (2). In the model (2), only the effects of covariates (W N y t , y t−1 , X t ) are allowed depending on τ. Obviously, when T is relatively small compared with N, it is difficult to estimate the distribution of individual fixed effects that depend on τ. Therefore, we restrict the estimates of η 0 to be independent of τ across all quantiles. Similar to Koenker [34] and Galvao [35], we treat the individual fixed effects as pure location-shift parameters to all quantiles and may be subject to shrinkage toward a common value as in the Gaussian random effects paradigm.
Koenker [34] proposed a general approach to estimating the quantile regression models for panel data. Using this method to the model (3), for fixed τ, the quantile regression estimators of (λ(τ), γ(τ), η 0 , β(τ)) can be derived by (λ,γ,η,β) = arg min (λ,γ,η,β) where ρ τ (u) = u(τ − 1{u ≤ 0}) is the check function, 1{·} is an indicator function, and v τ is a nonnegative scale weight used to control the influence of different quantiles on the estimation of η 0 ,ȳ it = ∑ N j=1 w ij y jt and η 0i = η 0 e i , where e i is an N × 1 vector in which the ith element is equal to 1 and the other elements are equal to 0.
In general, the estimators that we obtained from the above method are usually biased due to the endogeneity of the model. This problem can be improved by using instrumental variables, which affectȳ it and y i,t−1 in an appropriate way and are independent of ε it . According to Chernozhukov and Hansen [36] and Galvao [35], suppose that q × 1 vector z it is an available instrumental variable, and for fixed τ, define the objective function Here, v τ is a positive scalar weight, which is not important in parameter quantile regression, and we can set it as 1 in practice. For a given quantile τ, we define the IVQR estimators for the model (3) as where η (λ, γ, τ),β (λ, γ, τ),δ (λ, γ, τ) = arg min with x A = √ x Ax, and A is a positive define matrix. The final parameter estimators of interest parameters are given by Intuitively, a valid instrumental variable z it , which is independent of ε it , should have a zero coefficient. Hence, the finite sample objective function Q NTτ (λ, γ, η, β, δ) meets the identification condition, which is that the estimate of δ is close to zero when (λ, γ, β ) is close to the true population values λ(τ), γ(τ), β (τ) . We look for parameter values for λ (τ),γ(τ),β (τ) through the inverse step (7) so that the value of coefficient δ(λ, γ, τ) on the instrument variable is driven as close to zero as possible. A grid research can be used to implement the IVQR procedure in practice: (i) For a given τ, define a pair of grid values {(λ j , γ l ) : |λ j | < 1, |γ l | < 1, j = 1, · · · , J, l = 1, · · · , L} and run the ordinary quantile regression of (y it − λ jȳit − γ l y i,t−1 ) on (e i , x it , z it ) to obtain η (λ j , γ l , τ),β (λ j , γ l , τ),δ (λ j , γ l , τ) = arg min (η,β,δ) Q NTτ (λ j , γ l , η, β, δ).

Remark 1.
The pair of grid values of (λ, γ) can be determined based on actual requirements: the smaller the grid values are divided, the more accurate the estimators can be, and of course, the more time it takes.
However, for this case, the amount of calculation depends on the number of estimated quantiles and is very large. Therefore, we can use a numerical optimization function in R instead of the grid research according to Galvao [35].

Asymptotic Properties
Let Λ, Γ, E , and B be the parameter spaces of λ, γ, η, and β, respectively. In order to prove the asymptotic properties of the IVQR estimators defined in the above section, we need to establish the following assumptions. Assumption 1. The disturbance terms {ε it } are independent of each other, covariance stationary, with conditional distribution functions F it , and differentiable conditional densities f it , 0 < f it < ∞ with bounded continuous derivatives f (1) it , where i = 1, · · · , N and t = 1, · · · , T. Assumption 2.
and column sums; (iii) The sequence of matrices {H −1 (λ, γ)} are uniformly bounded in either row or column sums, uniformly in (λ, γ) ∈ Λ × Γ; (iv) The elements w ij of spatial weights matrix W N are uniformly at most of order l −1 N , such that l N → ∞ and l N /N → 0 as N tends to infinity.
Remark 3. Assumption 1 limits the density function of ε it . Assumption 2 is some essential features of the spatial weights matrix in the literature. Assumption 2 (i) guarantees that the random disturbances are well defined. By Kelejian and Prucha [40], and Lee [41], Assumption 2 (ii) limits the spatial correlation to a certain degree but helps to derive the asymptotic properties of spatial parameter estimator. Similar to Su and Yang [42], Assumption 2 (iii) requires {H −1 (λ, γ)} to be true uniformly in (λ, γ) ∈ Λ × Γ. Assumption 2 (iv) requires that w ij tends to 0 uniformly as N → ∞, which is also adapted from Su and Yang [42]. Assumption 3 emphasizes the global identifiability. Assumption 3 (i) requires the influence of instrument Z on the conditional distribution of y at relevant points, similar to Chernozhukov and Hansen [36]. Assumption 3 (ii) imposes the parameter space is compact and convex. Assumption 3 (iii) requires that the image of Λ × Γ × E × B under the mapping (λ, γ, η, β) → R τ (λ, γ, η, β) is simply connected. By Chernozhukov and Hansen [43], and Galvao [35], the mapping is a homeomorphism between (Λ × Γ × E × B) and R τ (Λ, Γ, E , B). Assumption 4 gives the condition of matrices for guaranteeing asymptotic normality. Assumption 5 imposes a restriction on the compactness of the parameter space of λ(τ), γ(τ) and requires that the objective function is convex at (λ, γ). Galvao [35] also assumed Assumption 6, which limits the boundaries of the variables. Assumption 7 is the same assumption as A6 in Galvao [35].
. Next, we study the consistency and asymptotic properties of the IVQR estimators defined in (7)- (9). All proofs are given in Appendix A. where

Remark 4.
The asymptotic covariance matrix Ω depends on the weights matrix A. If dim(δ) = dim(λ, γ), then the choice of A does not affect the result, and we can simply choose it as the identity matrix. According to Galvao [35], if dim(δ) > dim(λ, γ), then the choice of A affects the efficiency, and we can choose the inverse of covariance matrix ofδ(λ, γ, τ).

Monte Carlo Simulations
In this section, Monte Carlo simulations are used to assess the finite sample performance of proposed estimators. We use the bias, mean square error (MSE), relative bias (rB), and relative mean square error (rMSE) as the evaluation criteria. Here rB = where mcn is the total number of simulations, ξ j is the jth (j = 1, 2, . . . , mcn) simulated estimate, and ξ 0 is the true value of parametric estimates. The results to be reported in this section were estimated using R routines.
We run a small simulation experiment with mcn = 5000 and use data generated from the following model: . The spatial weightsmatrix W N is generated under Rook contiguity and Case group interaction, and a detailed description can be found in Su and Yang [42] and in Lin and Lee [45]. We consider N = {30, 50, 80}, T = {5, 10, 15}, and the values for remaining parameters as follows: We always discard the first 50 observations of y it to ensure that x it is not affected by its initial values.
In order to investigate the Monte Carlo results with variations on different quantiles, we try τ = 0.25, 0.5, and 0.75. According to Su and Yang [42], and Galvao [35], we set the instrumental variables as ∑ j w ij x jt and x i,t−1 . Our simulation results are presented in Tables 1-6 under three cases By observing the simulation results in Tables 1-6, we find the following conclusions: (1) The biases, MSEs, rBs (in percent), and rMSEs (in percent) are all small whether the random error distribution follows symmetrical distributions (N(0, 1) and t (3)) or a asymmetric distribution (χ 2 (3)). The MSEs and rMSEs of estimates drop quickly with increasing T. The accuracy of estimates improves significantly with increasing N or T. (2) The results in Tables 1 and 3 show that, under normal and t distributions of random error term, the rMSEs of estimates at quantile τ = 0.5 are smaller than that at quantiles τ = 0.25 and τ = 0.75. Tables 2 and 4 display the same phenomenon. The results in Table 5 show that, under χ 2 distribution of the random error term, the rMSEs of estimates at quantile τ = 0.25 are the smallest, followed by the quantile τ = 0.5, and finally, the quantile τ = 0.75. Table 6 displays the same phenomenon.         Tables 1 and 3, it is found that the rMSEs of estimates under normal errors are smaller than that under t errors, what implies that the accuracy of estimates under normal distribution of random error term is higher than that under t distribution; Tables 2 and 4 display the same phenomenon. Comparing Tables 1 and 5, it is found that the rMSEs of estimates under normal errors are smaller than that under χ 2 errors except for the special case of quantile τ = 0.25 and T = 5; Tables 2 and 6 display the same phenomenon. Comparing Tables 3 and 5, it is found that the rMSEs of estimates under t errors are smaller than that under χ 2 errors expect for lower quantile; Tables 4 and 6 display the same phenomenon. (4) Under the same distribution of random error term, the accuracy of the estimates has no significant difference for Rook and Case spatial weights matrices. This indicates that our estimation method is less affected by the choice of the spatial weights matrix.

Application
The tourism industry plays an important role for promoting economic development, social employment, and culture exchanges, and it has been paid much attention by governments in the world. As our model studies not only the influencing effects of covariates but also the spatial effects and dynamic effects of response variable, the lagged response variable, covariates, and the lagged covariates, it is interesting to use its estimation method to explore the influencing factors of international tourism foreign exchange earnings in China.
In this application, we mainly consider the influencing factors of international tourism foreign exchange earnings from the aspects of living standard, trade openness, tourism resources, transportation convenience, and service facilities in China. The data are based on a panel data of 31 provinces in China from 2011 to 2017 and were collected from the China Statistical Yearbook (http://www.stats.gov.cn/tjsj/ndsj/, 24 October 2018) and The Yearbook of China Tourism Statistics (https://data.cnki.net/yearbook/Single/N2020030 028, 1 November 2018). The established model is given by where response variable y t is provincial international tourism foreign exchange earnings and x t = log(PGDP t ), log(TO t ), log(TA t ), log(TC t ), log(SH t ) . Covariates PGDP t , TO t , TA t , TC t , and SH t are measured by the following methods: PGDP t : per capita GDP. It represents the living standard of provincial residents. TO t : proportion of the total value of imports and exports trade in GDP multiplied by 100. It is used to describe the provincial trade openness.
TA t : total number of A-rated tourist attractions. It reflects the level of provincial tourism resources.
TC t : weighted panel standardized indicators of civil aviation throughput, highway network density, railway network density, and inland waterway network density. According to the suggestion of Zhou and Sheng [46], we obtained the weighted panel standardization index of above four indicators by following formula where z 1,it , z 2,it , z 3,it , and z 4,it denote the observations of civil aviation throughput, highway network density, railway network density, and inland waterway network density, respectively; i = 1, · · · , 31; t = 1, · · · , 7; and k = 1, · · · , 4. According to the importance of these four indicators, we assigned weights ω 1 = 0.3, ω 2 = 0.3, ω 3 = 0.2, and ω 4 = 0.2. TC t is used to reflect provincial transportation convenience. SH t : total number of star-rated hotels. It represents the level of provincial service facilities.
The spatial weights matrix we adopted is the inverse of the greater circle distance calculated based on the longitude and latitude of China's provincial capital cities. It is standardized to ensure that the sum of each row elements is equal to 1. The estimation results of model (10) at quantiles τ = 0.10, 0.25, 0.50, 0.75, and 0.90 are reported in Table 7. For a comparative analysis, we also report the ordinary least squares (OLS) estimates in the last column of the Table 7.
According to Table 7, we find the following results: (1) For IVQR, all estimates of parameters are significant at test level 0.01 and not the same at different quantiles.
(2) For IVQR, the serial dynamic effects of response variable and covariates exist between adjacent years. The spatial spillover effects of response variable, covariates, and their lagged terms exist among Chinese provinces. (3) The estimates of coefficients for covariates log(PGDP t ), log(TO t ), and log(TC t ) are uniformly positive at different quantiles, which implies that these influencing factors have positive effects on international tourism foreign exchange earnings. The estimates of coefficient for covariate log(TA t−1 ) is uniformly negative at different quantiles, which indicates that this influencing factor has an inhibitory effect on international tourism foreign exchange earnings. (4) The estimates of the coefficient for the Durbin term W N log(TC t ) is uniformly positive at different quantiles, which means that the improvement of transportation convenience in neighboring regions has a positive impact on the provincial international tourism foreign exchange earnings. The estimates of the coefficients for the Durbin terms W N log(TA t ), W N log(TO t−1 ), and W N log(SH t−1 ) are uniformly negative at different quantiles, which show that the improvement of tourist attractions, trade openness, and star-rated hotels in neighboring regions all have negative effects on response variable. (5) The influences of other covariates on international tourism foreign exchange earnings have positive and negative effects at different earning levels. The estimates of most covariates at quantiles τ = 0.25, 0.50, and 0.75 have same positive or negative influence directions. (6) OLS estimates reflect the influence of conditioning covariates on the average of response variable, and they are quite different from the IVQR estimates. Variables log(y t−1 ), log(PGDP t ), log(TO t ), log(TC t ), W N log(TC t ) and W N log(TA t−1 ) have significant positive effects on international tourism foreign exchange earnings; variables W N log(PGDP t ), W N log(TO t ), W N log(TA t ), log(TO t−1 ), log(TA t−1 ), and W N log(TO t−1 ) have significant negative effects on international tourism foreign exchange earnings; and the influences of variables W N log(y t ), W N log(y t−1 ), log(TA t ), log(SH t ), W N log(SH t ), log(PGDP t−1 ), log(TC t−1 ), log(SH t−1 ), W N log(PGDP t−1 ), W N log(TC t−1 ), and W N log(SH t−1 ) on international tourism foreign exchange earnings are not significant. By establishing a quantile regression SDDPD model with fixed effects, we explored the influencing factors of international tourism foreign exchange earnings at different quantiles. In order to increase international tourism foreign exchange earnings, local governments should take appropriate measures to exert positive effects and to restrain the negative effects of influencing factors according to their development stage of foreign tourism. Variables log(PGDP t ), log(TO t ), and log(TC t ) have positive effects on international tourism foreign exchange earnings at different quantiles, which implies that all provinces should improve living standard, promote trade openness, and provide convenient transportation. Provinces at quantile 0.1 should appropriately control the number of tourist attractions and improve their quality. For provinces at other quantiles, the number of tourist attractions can be increased. Provinces at quantiles higher than 0.75 should appropriately control the number of hotels. For provinces at other quantiles, the number of hotels can be increased. Neighboring provinces should carry out coordination and cooperation, share tourism resources and facilities, avoid vicious competition, and jointly promote increases in international tourism foreign exchange earnings.

Conclusions
This paper focuses on studying estimation and inference in a quantile regression SDDPD model with fixed effects. Conventional fixed effects estimators of quantile regression specification are usually biased in their presentation of lagged response variable in spatial and time as regressors. We first simplified this model as a spatial autoregressive dynamic panel data model with fixed effects. In order to eliminate the endogeneity, we used the IVQR method proposed by Chernozhukov and Hansen [36] to obtain the quantile estimators of parameters. Under some regular assumptions, the consistency and asymptotic properties of IVQR estimators were derived. Monte Carlo simulations show that our estimators not only have outstanding finite sample performance but also have good robustness for different spatial weights matrices and for different disturbance term distributions. Finally, we illustrated our method with an application to analyze the factors affecting international tourism foreign exchange earnings of 31 provinces in China from 2011 to 2017, and the results we obtained are consistent with the actual situations.
Lemma A2. Under Assumptions 1-7, we have Proof of Lemma A2. Let π = (λ, γ, η, β, δ) be an arbitrary value in parameter space The relationship betweenε it (τ) and ε * it (τ) can be expressed as We adopt the idea of Galvao [35] to prove this lemma, with the major difference lying in the appearance of spatial dependence in our model. Consider partitioning the parameter space Π into finite disjoint subsets Π 1 , · · · , Π K N such that the diameter of each subset does not exceed r NT = N a 4K N CT , where C is a constant. Let π j = (π j,η , π j,β ) be the fixed points in Π j , j = 1, · · · , K N . Then, the left-hand side of (A1) does not exceed U 1 + U 2 , where r,s . Let ς 1 = b r ε it −h r,s ε jl and ς 2 = b s ε jl −h s,r ε it . It is easy to obtain that the joint density of (ς 1 , ς 2 ) conditional on F is given by 1 |B r,s | f it b s ς 1 +h r,s ς 2 B r,s f jl b r ς 2 +h s,r ς 1 B r,s , where B r,s = b r b s −h r,shs,r . Obviously, B r,s = 0 which is satisfied for sufficiently large N because according to Assumption 2 (ii) and two evident facts of Kelejian and Prucha [40], h r,shs,r = O(l −2 N ) = o(1) and b r b s → 1. In addition, the marginal density of ς 1 is given by 1 |b r | f it ς 1 +h r,s ς 2 b r f jl (ς 2 )dς 2 .
Proof of Theorem 2. Consider a set of closed balls B n λ(τ), γ(τ) with radius r n centered on λ * (τ), γ * (τ) , where r n → 0 slowly enough. Define wherex it = (ȳ it , y i,t−1 , x it , z it ) as defined in Lemma A2. For a fixed d nβ , we first consider the behavior of d nη e i for each i. Define For a fixed d nβ and d nδ , λ n (τ) − λ * (τ)