{'fraction_non_alphanumeric': 0.03512655182369358, 'fraction_numerical': 0.01818881724569339, 'mean_word_length': 4.586827956989247, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Nonparametric Instrumental Variable Regression Approach to Quantile Inference

This study extends the Bayesian nonparametric instrumental variable regression model to determine the structural effects of covariates on the conditional quantile of the response variable. The error distribution is nonparametrically modelled using a Dirichlet mixture of bivariate normal distributions. The mean functions include the smooth effects of the covariates represented using the spline functions in an additive manner. The conditional variance of the second-stage error is also modelled using the spline functions such that it varies smoothly with the covariates. Accordingly, the proposed model allows for considerable flexibility in the shape of the quantile function while correcting for an endogeneity effect. The posterior inference for the proposed model is based on the Markov chain Monte Carlo method that requires no Metropolis-Hastings update. The approach is demonstrated using simulated and real data on the death rate in Japan during the inter-war period.


Introduction
A set of quantiles provides a more complete description of a distribution than the mean. Similarly, exploring the relationship between the quantiles of the conditional distribution of the response variable and a set of covariates offers key insight into the structure of the data at hand. As an important approach to quantile inference, the quantile regression has received substantial attention since Koenker and Bassett's (1978) seminal work. There is a vast literature on the theory of quantile regression; for an overview, see Koenker (2005), Yu et al. (2003), and Buchinsky (1998).
Studies on Bayesian quantile regression particularly increased following Yu  As in the case of conventional mean regression, a standard quantile regression estimator is known to be biased when an endogenous variable is included. The problem of endogeneity in the quantile regression framework was recognised as early as in the 1980s (Amemiya, 1982;Powell, 1983).
Since then, there has been a growing interest in the inference for quantile regression models with endogenous variables and quantile treatment effects in the context of the frequentist approach (e.g. However, despite the growing interest in and demand for Bayesian quantile inference, the literature on the Bayesian approach to quantile inference in the presence of an endogenous variable remains particularly sparse. Lancaster and Jun (2010) adopted the moment conditions adopted in Chernozhukov and Hansen (2006) for the empirical likelihood. Drawing on Lee (2007), Kobayashi (2016) and Ogasawara and Kobayashi (2015) employed the control function approach. They considered the use of parametric and semiparametric first-stage error distributions, whose quantile is restricted to zero for some quantile level, and employed the asymmetric Laplace distribution in the second stage, as in the standard Bayesian quantile regression. Kobayashi's (2016) simulation study showed that while this approach works well in many situations, it may provide biased estimates owing to the restrictive specification of the first-stage error distribution. Furthermore, both approaches assume a simple linear model for the conditional quantile and monotonicity (non-crossing) of the estimated quantiles is not necessarily guaranteed.
Given the preceding discussion, this study aims to provide a flexible Bayesian approach to determine the structural effects of covariates on the conditional quantile of the response variable in the presence of an endogenous variable. To this end, we extend the Bayesian nonparametric instrumental variable regression model proposed by Wiesenfarth et al. (2014) and Conley et al. (2008) and consider a location-scale specification. In particular, the joint error distribution of the instrumental variable regression model is modelled using the Dirichlet process mixture of bivariate normals to capture the non-normality of the error distribution. To capture the nonlinear effect of the covariate on the response variable, the mean functions are nonparametrically modelled using the spline functions.
In addition, we also nonparametrically model the conditional variance of the second-stage regression using spline functions such that the conditional variance varies smoothly with covariates. Unlike the quantile curves implied from the Bayesian nonparametric instrumental variable regression model, those implied from the proposed model are not necessarily parallel to each other since the conditional variance differs by covariate value. Therefore, the estimated quantile curves allow significant flexibility in estimating the structural effect of the covariate of interest. In addition, since the conditional quantiles are computed from the nonparametrically estimated error distribution, the monotonicity of the quantile curves is guaranteed. The proposed model is estimated using the MCMC method that includes no Metropolis-Hastings (MH) update. The MCMC method is applied to sample from the posterior distributions of the conditional quantiles.
The remainder of this paper is organised as follows. The proposed model is introduced and the choice of the prior distributions is discussed in Section 2. The MCMC method for the posterior inference on the conditional quantiles is developed in Section 3. The proposed approach is illustrated using the simulated data in Section 4 and real data in Section 5, where the effect of doctors on the death rate in Japan during the inter-war period is studied. Finally, we conclude in Section 6.

Bayesian Instrumental Variable Regression Model
The following parametric Bayesian instrumental variable regression model provides the starting point of the model proposed in this study: where y is the response variable, d is the endogenous variable, z is the instrumental variable, γ 0 and β 0 are the intercepts, γ 1 and δ are the coefficient parameters, and = ( 1 , 2 ) is the error term following The correlation coefficient ρ = σ 12 / √ σ 11 σ 22 represents the degree of endogeneity.
The normality assumption for can be restrictive. Outliers, some form of model misspecification, and heterogeneity can also induce non-normality in the error term. To relax the normality assumption, where γ 0 and β 0 are intercepts and f (z) and g(d) represent the nonlinear effects of the instrument and endogenous variables, respectively. The assumptions required for identification are denoted by The smooth functions are approximated using the spline functions and are represented as the linear combinations of the basis functions. The extensive simulation in these studies revealed that the nonparametric modelling approach is more robust against non-normality and outliers and performs better than existing approaches, including two-and threestage least squares, a model under the normality assumption, and the control function approach with the generalised cross-validation procedure for the choice of a smoothing parameter.

Nonparametric Model for Quantile Inference
To allow for inference on the conditional quantile while correcting for an endogeneity effect, we extend (2) and consider a location-scale model for more flexibility in the quantile curve. Furthermore, to introduce the model for the conditional variance and prior distribution in a more flexible manner, the error term is rewritten following Lopes and Polson (2014). Specifically, this study considers the model given by for i = 1, . . . , n, where α 0 is the intercept and h is the smooth function of d. For identifiability, the functions f (·), g(·), and h(·) are centred around zero. This formulation implies that the components of the covariance matrix Σ i for ( 1i , 2i ) is recovered through σ 11i = τ i , σ 12i = η i σ 11 i , and σ 22i = s i σ i + σ 2 12i /σ 11i . Although each equation in (3) includes only one nonlinear effect for notational simplicity, the model can include more functions of covariates in an additive manner.
As in Conley et al. (2008) and Wiesenfarth et al. (2014), the error distribution is nonparametrically modelled using the Dirichlet process mixture, but the prior distribution is introduced in a different manner. We introduce the Dirichlet process prior for the distribution of θ i = (µ 1i , µ 2i , η i , τ i , σ i ): where DP(a, G 0 ) denotes the Dirichlet process with the precision parameter a and base measure G 0 .
Since our model contains the global intercepts γ 0 and β 0 , to achieve identifiability, we employ the priors µ 1 and µ 2 such that µ 1 and µ 2 have mean zero apriori and force n i=1 µ 1i = n i=1 µ 2i = 0 to ensure they have zero mean posteriori. In addition, since s i includes the intercept α 0 , we force  (Lopes and Polson, 2014). In this study, we employ the following independent distributions for G 0 : τ ∼ IG(t 1 , t 2 ), and σ ∼ IG(s 1 , s 2 ). The choice of the hyperparameters and hyperpriors is discussed in Section 2.3.
Each smooth function of covariate is approximated using a linear combination of B-spline functions in the study. For f (·), for example, given the number of interior knots K f and degree m f , let us denote the set of knots by can be calculated through a simple recursion of de Boor (2000). The smooth function f (z) is approximated by Similarly, functions g and h are approximated using b Kg,mg (d)β 1 and b K h ,m h (d)α 1 . The smoothness of the functions is controlled through the prior distributions on the basis coefficients by penalising the differences between adjacent coefficients (see Section 2.3).
Finally, we are interested in estimating the structural effect on the conditional quantile: where Q e|d (p) is the quantile function of e. In addition, since the quantiles of the error term are directly obtained from the nonparametrically modelled distribution, crossing of estimated quantiles is not an the conditional mean. This is different from the Bayesian quantile regression approach of Kobayashi (2016) and Ogasawara and Kobayashi (2015), who employed Lee's (2007) control function approach in the quantile regression framework. The control variable is introduced to the second-stage regression such that the p-th quantile for some p of interest is corrected to be equal to zero. Then, assuming a linear model, the structural effect is estimated fromS y,d (p) = β 0 + dδ. While this approach works well in several situations, the required assumption for the first-stage error that the α-th quantile is equal to zero for some α ∈ (0, 1) can be quite restrictive in modelling the error distribution. In terms of the quantile modelling, the present approach may be considered a Bayesian nonparametric extension of Gilchrist (2008), who considered various ways to introduce the covariate influence on the quantile function using parametric models.

Default Prior Distributions
This section discusses the choice of the prior distributions and associated hyperparameters. For the spline coefficients for the nonlinear effects-γ 1 , β 1 , and α 1 -the Gaussian random walk priors are assigned following Lang and Brezger (2004). In the case of γ, the joint prior density is given by For the intercepts γ 0 , β 0 , and α 0 , we adopt the normal priors N (0, V γ 0 ), N (0, V β 0 ), and N (0, V α 0 ), respectively. Since we do not have information on these parameters, the hyperparameters are set such that the priors are vague: We turn to the prior specification for the components in the Dirichlet process. Specifically, it is natural to assume that ϑ e is centred around this value. Therefore, ϑ e ∼ N (e 0 , T e0 ) is assumed, where e 0 is equal to the sample covariance between d and y divided by the sample variance of d and T e0 = 10. Since we do not know the variability of η, that is, the degree of heterogeneity in the endogeneity, we place a prior distribution for ϕ e . Namely, ϕ e ∼ IG(e 1 , e 2 ) with e 1 = e 2 = 2 such that ϕ e takes a value between 0 and 3 with high probability. For τ and σ, assuming that the data are appropriately rescaled, we set t 1 = t 2 = 2, s 1 = 2, and s 2 = 1. Finally, the precision parameter of the Dirichlet process follows the gamma prior a ∼ G(a 1 , a 2 ) with a 1 = a 2 = 2, such that both small and large numbers of mixture components are allowed. In the simulation study in Section 4, we also consider some alternative prior specifications for ϑ e and a.
with ω k ∼ B(1, a), and B(a, b) is the beta distribution with parameters a and b (Sethuraman, 1994).
We also let k * denote the minimum integer such that k * k=1 > 1 − min {u 1 , . . . , u n }. To apply the slice sampler and retrospective sampler to the variables involved in the Dirichlet process, we work on the following joint density: where N (·; µ, σ) denotes the density function of the normal distribution with mean µ and variance σ, In addition, collecting terms of (3) and (4) would be useful for some steps of the Gibbs sampler: Our Gibbs sampler proceeds by alternately sampling from the full conditional distributions of : Generate k i from the multinomial distribution with probabilities • Sampling a: By introducing r ∼ B(a + 1, n), the full conditional distribution of a is a mixture of two gamma distributions given by where n * is the number of distinct clusters and ϕ/(1 − ϕ) = (a 1 + n * − 1)/(n(a 2 − log r)). See Escober and West (1995).
• Sampling ϑ 1 and ϑ 2 : The full conditional distribution of ϑ 1 is N (θ 1 ,T 1 ), wherê The sampling for ϑ 2 is performed in a similar manner.
• Sampling ϑ e and ϕ e : The full conditional distribution of ϑ e is N (θ e ,T e ), wherê and the full conditional distribution of ϕ e is IG(ê 1 ,ê 2 ), wherê • Sampling α: Note that where e * i = log e 2 i follows the log-χ 2 distribution with one degree of freedom. It is known that the log-χ 2 distribution can be accurately approximated by the ten-component mixture of normal and then, sampling α from N (α,Â), wherê The values for ξ h , ζ h , and h are presented in Appendix A. See also Chib and Greenberg (2013).
• Sampling φ f , φ g , and φ h : We sample φ g and φ h in the same manner.

Quantile Estimation
The t-th iteration of the Gibbs sampler provides a posterior draw from G: To estimate (5) for any p ∈ (0, 1), we obtain draws x β (t) +Q where Φ(·) denotes the distribution function of the standard normal distribution. As in Taddy and Kottas (2010), F (e|d) and Q e|d (p) must be computed over a sufficiently fine grid of covariate values.
For a high dimensional covariate vector, the inference would be feasible only when a one-or twodimensional grid is used for the covariates of interest, while the other variables are fixed at some value.
This specification for the mixture components implies that the correlation coefficients between v and e are 0.894, 0.287, and 0.624. The data are replicated 100 times.
We fit a simple version of the proposed model comprising a linear mean and variance functions but with a nonparametric error distribution. For comparison, we also fit two alternative models.
posterior mean of (5) at d g for the r-th replication. We set G = 100 and let d 1 = max r {min i {d i }}. Figure 2 shows that the biases for the restricted and uncorrected models nonlinearly depend on the value of d and do not become smaller as the sample size increases.
By contrast, for all quantiles, the biases for the proposed model are near zero for the region considered and become closer to zero as the sample size increases. Similar to the biases, Figure 3 shows the clear dependence of RMSEs for the restricted and uncorrected models on the value of d.

Nonparametric Models
In this section, the proposed model is demonstrated under various settings for the mean and variance functions. The data are generated from where we use the same setting for the error distribution as in Section 4.1. For the mean and variance functions, the following three settings are considered: (i) f (z) = 1.5z + 2 exp(−16z 2 ), g(d) = 1 + sin(d), s(d) = exp(0.5 − 0.3d 2 ),   As shown in Figure 6, the RMSE for the proposed model is smaller than those for the restricted and uncorrected models for most parts of the region where d is observed. Here as well, the RMSEs for the restricted model and, particularly, the uncorrected model, substantially vary by the value of d.
The figures also show that there is little difference in the biases and RMSEs between the default prior specification and with the two alternative prior specifications.  as the sample size increases. We obtained similar results for the RMSE, which are also shown in Figure 6. Furthermore, the results for these settings under the proposed model do not appear to be significantly influenced by the prior specification.

Background
To illustrate the proposed approach using real data, the effect of the number of doctors on Japan's death rate during the inter-war period is studied. The economic history literature has extensively studied declines in the mortality rates of Western countries across the 19th and 20th centuries. Recent studies have found the improving effects of national health insurance systems on mortality rate (e.g. Winegarden and Murray, 1998;Bowblis, 2010). The successful application of germ theory has also been considered to decrease the mortality rate (Mackenbach, 1996;Mokyr and Stein, 1997).
To the effect of medical care, the importance of opportunity for medical care, for instance, through one-on-one contact and medical vouchers, for people with limited access to medical care have also been stressed in recent studies (Moehling and Thomasson, 2014;Ogasawara and Kobayashi, 2015).
These studies revealed that increased access to medical care contributed to the decline in mortality rates in the process of economic development.

Data and Model
To analyse the effects of an increase in medical opportunities on death rate, we compiled a large set  doctors.
The covariates include the share of infant (child) and elder (old) population and the age structure of the population is assumed to have nonlinear effects on the death rate. We also include the log odds ratios of the proportion of population employed in the agricultural (agri), mining (mine), and industrial (ind) sectors and the log odds ratio of the primary school enrolment rate (enrol). These variables are included to control for the population and industrial structure effects, potential wealth level, and unobservable macroeconomic shock in a given year and are assumed to have linear effects.
To summarise, the following model is fitted to the data: where x l = (agri, ind, mine, enrol) and D(t), t ∈ {1925, 1930} is the dummy variable such that iterations. Every tenth Gibbs draw is retained for posterior inference.

Results
First, the results for error density and mean and variance functions are presented. Figure 9 presents the estimate for the joint error density. The figure shows that the error density exhibits a slight deviation from normality. Figure 10 presents the posterior means and 95% credible intervals for the nonlinear effects in the mean and variance functions. It is shown that lnarea is able to explain the variation in lndoc and appears to be a valid instrument. Moreover, the estimate for f 3 exhibits some nonlinearity.
From the estimate of g 3 , our key variable lndoc has some nonlinearly improving effect on the mean death rate. The estimate of g 4 shows that the mean death rate for the counties with a large number of doctors may be lower in 1930 than in 1925, although the 95% credible interval includes zero for the entire range of lndoc. The conditional variance of the death rate appears to vary nonlinearly with lndoc as well. The estimate of h 3 shows that the conditional distribution of the death rate has a large variance for the small values of lndoc and decreases as lndoc increases. Table 2 presents the posterior means, 95% credible intervals, and inefficiency factors for the linear effects.    can be attributed to the fact that our data include few observations in this region during 1930 (see Figure 11). Therefore, we limit our focus to region lndoc > 2.
As for the effect of doctors, the figure clearly shows that the quantiles for p = 0.5, 0.9, and 0.99 We also estimate the joint quantile surfaces over child and lndoc denoted by S dr|child,lndoc,xt,D(t) (p) and over old and lndoc denoted by S dr|old,lndoc,xt,D(t) (p).

Conclusion
In this study, we considered a Bayesian nonparametric approach to inference for the conditional

Acknowledgements
This study was supported by JSPS KAKENHI [Grant Numbers 15K17036 and 16K17153]. The computational results were obtained using Ox version 6.21 (Doornik, 2007).

Appendix A Parameters for Mixture Sampler
In the mixture sampler to sample α, the log-χ 2 distribution is approximated using the ten-component mixture of normal distributions: f (e * ) ≈