# The dynamics of the prices of the companies of the STOXX Europe 600 Index
through the logit model and neural network

Federico Mecchia  Master in Economics and Finance, Italy  Marcellino Gaudenzi
Dipartimento di Scienze economiche e statistiche, Università di Udine, via
Tomadini, 30/A, Italy

###### Abstract

The aim of the present work is analysing and understanding the dynamics of the
prices of companies, depending on whether they are included or excluded from
the STOXX Europe 600 Index. For this reason, data regarding the companies of
the Index in question was collected and analysed also through the use of logit
models and neural networks in order to find the independent variables that
affect the changes in prices and thus determine the dynamics over time.

Keywords: STOXX600, Logit model, Neural network, Stock prices

##  1 Introduction

Financial indexes today play the role of protagonists in the worldwide economy
and, in fact, they are often taken into consideration to acquire an immediate
and well-rounded sight about a specific market or about particular trends. In
addition to this, because of the important impact of these indexes on the
global economy, determining the dynamics of the prices of the stocks of the
companies of the various indexes can be of fundamental importance, also with
regard to the assessment of other correlated dynamics and the impact on the
economic world as a whole.

For these reasons, the focus of this work is understanding the dynamics of the
prices of stocks of the companies of the STOXX Europe 600 Index; for this
purpose, the logit model and the neural network have been taken into
consideration and applied, through the ad hoc code written in Python, to the
created database in order to reach the aforementioned objective.

As it was stated earlier, indexes play a key role in today’s economy and
various researches have been carried out with regard to this topic. For
example, researchers have investigated different aspects of the S&P500 Index,
of the NASDAQ Index, of the DOWJONES Index and of the NIKKEI Index. To this
regard, in fact, for example, Lee and Hao [13] have considered the S&P500
Index in order to study the asymmetric effects of error corrections between
oil prices in the U.S.A. and the prices of the aforementioned index, while
Lento and Gradojevic [14] focused their attention on understanding the price
spillovers of the S&P500 Index related to the effects on the market of
COVID-19. Moreover, Metghalchi et al. [16] have explored the trading
strategies related to the moving average with regard to the NASDAQ Index,
Giannarakis et al. [9] have carried out an analysis of the Dow Jones
Sustainability Index and Christiansen and Koldertsova [7] have taken into
account several indexes, such as the NASDAQ Index, in order to explore the
role of stock exchanges in corporate governance.

With regard to the STOXX Europe 600 Index, Belhouchet and Chouaibi [2] took
into consideration the STOXX Europe 600 Index to determine the relationship
between the integrated reporting quality and the audit committee attributes
and Kyriazis [10] analysed the link between COVID-19 deaths, crude oil and
gold and the STOXX Europe 600 Index and its sub-indices and STOXX Europe 50.
In addition to this, Salvi et al. [19] studied the relationship between the
volume of LBOs and the STOXX Europe 600 Index and De Blasis and Petroni [8]
considered the Index in question to analyse the relationship between oil and
renewable energy firms with regard to the COVID-19 pandemic. In addition to
this, Langenstein et al. [12] took into consideration the companies of the
STOXX Europe 600 Index together with the companies of the S&P500 Index to
explore the best strategies that can be enhanced during crisis, Chen et al.
[6] analysed the dynamics of the FTSE 100, S&P500 Index and STOXX Europe 50 by
taking into account different exchange rates and Bonsón and Bednárová [3]
considered some of the companies of the STOXX Europe 600 Index to carry out
their research regarding the CSR reporting practices of the companies of the
Eurozone. To the authors’ knowledge, there are no researches focusing on
whether the inclusion or the exclusion of a company from the STOXX Europe 600
Index has an impact on the fluctuations of the prices of the company in
question; for this reason, the objective of this work is understanding the
dynamics of the prices of companies depending on whether the company
considered is a constituent or is not part of the STOXX Europe 600 Index.

For this purpose, several datasets were built by taking into consideration the
data that can be found in Bloomberg. These datasets were then analysed through
the logit model to determine which are the most significant variables that
have an impact on the fluctuations of the prices. The selected variables were
then taken into account to build the neural networks for each of the companies
taken into account; the levels of accuracy varied according to the company
considered, but overall it is possible to state that they were all around the
value of 70%, thus providing a sufficiently satisfactory result.

If compared to the existing literature, the present work is characterized by
some distinctive features. Among these, there is the fact that the logistic
regression model is used to select the independent variables that are then
included in the creation of the neural network and therefore the logit model
and the neural networks are not considered to be independent procedures, but
part of the whole process. Moreover, another distinguishing aspect of this
work is the thorough process of data selection and data analysis that have
been carried out. Lastly, another peculiar characteristic of this work is the
comparison of the neural networks among various companies of the STOXX Europe
600 Index in order to determine the overall accuracy derived from the
application of neural networks to different realities of the Index in
question.

##  2 Methodology

As it was stated earlier, the methodology employed in this work in order to
analyse the dynamics of the prices of the companies of the STOXX Europe 600
Index involved the use of the logit model and of the neural networks, that
have been applied to the ad hoc created database through the code written in
Python. More specifically, the logit model was employed in order to determine
which are the most significant independent variables that have an impact on
the fluctuations of prices of the companies considered; these variables were
then taken into account to build the neural networks for each of the companies
taken into consideration.

The logit model is a non-linear regression model and it is used in order to
determine the probability that a particular event occurs. The dependent
variable that is used in this type of model is dichotomic; in this case the
dependent variable is represented by the price fluctuation, as it is explained
in the following parts of this work. As it was stated earlier, understanding
the dynamics of the fluctuations of prices is of fundamental importance,
because prices are among the most important characterizing features of data
concerning stocks and also because the impact they can potentially have on
other aspects of the economy is a key element of various financial scenarios.

In technical terms, for each company taken into account consider the dependent
dichotomic variable $Y_{i}$, which therefore can only be equal to either 0 or
1 and which indicates whether there has been a negative or positive change in
price with respect to the prior period, with $i$ ($i=1,\dots,n)$ indicating
the specific time instant taken into consideration. Moreover, consider the
vector $x_{i}=(x_{ij},\dots,x_{im})$ and which includes at the time instant
$i$ all the types of data regarding the various characteristics of the company
taken into account, indicated by $j$ (with $j=1,\dots,m)$. The logit model can
therefore be defined as:

$logit(P(Y_{i}=1|x_{ij}))=\log\left(\frac{P(Y_{i}=1|x_{ij})}{P(Y_{i}=0|x_{ij})}\right)=\log\left(\frac{P(Y_{i}=1|x_{ij})}{1-P(Y_{i}=1|x_{ij})}\right)=\beta_{i}$

Furthermore, the neural networks can be thought of as the artificial
transposition of biological systems of neurons that contribute to the learning
mechanisms of humans. In other words, as the following image shows (Figure 1),
the neurons are structured according to an architecture so that each neuron is
connected to the other neurons of the network in question. In addition to
this, every neural network has different layers. The first layer is known as
the input layer, the last layer is known as the output layer and the layers in
the middle are the hidden layers.

Figure 1: Architecture of neural network

Each neuron of a specific layer performs different functions. First it sums up
all the values of all the neurons of the previous layer to which it is
connected and multiplies these values by the weight, which is located between
the previous neuron and the neuron in question; in addition to this, the bias
value can also be added. The neuron then uses the activation function which
converts the predetermined value to a number between 0 and 1. This process is
repeated throughout all the neural network, until the final layer (output
layer) is created. Finally, if the final output matches the expected output,
it means that the computed parameters are valid and so the process can be
carried out with regard to the other inputs; in all other cases the weights
are changed and the process is repeated until a satisfactory result is
achieved.

With regard to the use of logit models for prediction purposes, various
authors have carried out their researches by applying this type of model to
different situations. To this regard, for example, Celine et al. [5] employed
the logit model for the prediction of employability, while Borucka [4] worked
with the logistic regression model with regard to the transport services.

In addition to this, other authors have focused their attention on the
application of the logistic regression model to the economic and financial
sectors. To this regard, for example, Zaini et al. [23] used the logistic
regression to analyse the market movements and Ali et al. [1] employed the
logistic regression for stock prediction of the PSX. In addition to this,
Smita [21] used this type of model for the prediction of the performance of
S&P BSE30 Company and Mironiuc and Robu [17] used the logit model to predict
the performance of stocks on the emerging markets.

With regard to the use of neural networks for price prediction, Mehtab and Sen
[15] employed different tools, including the neural networks, for the
prediction of stock prices and Shahvaroughi Farahani and Razavi Hajiagha [20]
studied how neural networks and metaheuristic algorithms can be used for stock
prediction. Lastly, Song et al. [22] used different types of neural networks
to predicts stock prices, Moghaddam et al. [18] also employed the neural
networks to analyse the NASDAQ stock price and also Lachiheb and Gouider [11]
worked with neural network for the prices prediction of the TUNINDEX stocks.

##  3 Construction of datasets and analysis of data

As it was stated earlier, Python was employed in order to carry out all the
different parts of the analysis and model creation. All the data that was
taken into consideration for this work was taken from Bloomberg. More
specifically, two sources of data were considered in order to create the final
databases that were later involved in the parts inherent to the data analysis.

To this regard, in fact, the first source of data consists of the list of the
600 companies that compose the STOXX Europe 600 Index for each Friday from
Friday 4th January 2002 until Friday 31st December 2021. In case the data for
the specific Friday taken into consideration was not available, the previous
day before the Friday in question for which the data could be collected was
taken as reference for that week. This happened in few cases; for example, the
data for Friday 09th April 2004 was not available and therefore the data of
Thursday 08th April 2004 was considered.

The second source of data, instead, consisted of historical data of each of
the companies that were recorded at least once in the aforementioned 20-year
period.

In the beginning the first source of data, therefore the lists of the
components of the STOXX Europe 600 Index were considered in order to carry out
the first phases of the analysis. More specifically, the dynamics regarding
the presence or absence of the companies throughout the 20 years period was
investigated. To this regard, the number of times a specific company was a
constituent of the STOXX Europe 600 Index greatly varies according to the
specific company taken into account. In order words, some companies were
present just a limited number of times throughout the 20-year period and
therefore not for the entire timespan.

The range regarding the number of times a company was a constituent of the
STOXX Europe 600 Index was divided into five groups, more specifically each
representing one fifth of the range. Each class, as it can be easily inferred
from Figure 2, was composed of different number of companies.

Figure 2: Number of companies in each group

A total of 50 companies, 10 from each of the aforementioned class, were then
randomly selected in order to carry out the other parts of the analysis.

As it was previously stated, the second source of data consisted of historical
data of each of the 50 companies that were randomly selected. The historical
data consisted of financial data, data regarding the member of boards, and
other types of data. These files containing the historical data for each
company were then merged with the files containing the list of the
constituents of the STOXX Europe 600 Index in order to identify when (in order
words the specific dates), a particular company was present or absent with
regard to the Index in question. For this purpose, a new column was created in
each file of each company considered and the column in question was filled
through the merge procedure with values equal either to 0 or to 1 (the value 0
indicates that the company was not part of the Index at the date, while the
value 1 indicates that was present in the Index at that date).

In addition to this, also another column was added to each dataset. The column
in question indicated the dynamics of the prices, in order words whether the
price moved up or down with regard to the previous observation. To this
regard, the value 1 was assigned whenever there was a positive change in the
price and the value 0 was used to fill the column whenever there was a
decrease in the price considered.

Moreover, some potential independent variables were considered with a 7-days
lag; in other words the values of the previous week of some variables were
taken into account and this was done to understand how the data of the
previous Friday (or of the previous week) could have an impact on the
dependent variable seven days later.

Through the selection of specific companies, it was possible to observe that
data regarding the 20-year timespan in some cases was lacking and therefore,
depending on the entity of the missing data, some variables were removed. In
addition to this, the timespan for each company was designed ad hoc on the
basis of the data available; this fact inevitably caused the reduction of the
total timespan considered to shorter periods of time depending on the company
considered.

The following part of the analysis then consisted in normalizing the various
aforementioned datasets. To this regard, the following formula was used in
order to normalize each data of the datasets:

$x_{\hbox{normalized}}=\frac{x-\min{x}}{\max{x}-\min{x}}$

Therefore, first the minimum and maximum values of each column were computed
and then the minimum value was subtracted from the specific data taken into
consideration and the result was then divided by the range, defined as the
difference between maximum and minimum values. In this way, all the data of
the various datasets were converted in values ranging from 0 to 1. Moreover,
the missing data was replaced by the mean of the column considered (in other
words the mean of the available data for the specific variable taken into
account).

The different datasets were then taken into consideration with regard to the
logistic models. The logistic models were thus built taking into account the
datasets of the various companies with the column regarding the price increase
or decrease as dependent variable and all the other columns as independent
variables.

By building these models it was possible to identify which were the most
significant independent variables that could be taken into account for the
next parts of the analysis. In fact, through the logit models it was possible
to observe that among the most significant variables there were the variable
indicating whether a specific company was present or absent at a specific date
in the Index in question, the values regarding the total return (gross
dividends) with a 7-day lag, the average sentiment value of the news of the
day and the number of trades of that day. It was therefore possible to notice
that the fact that a company belongs to the Index in question at a specific
time instant has an impact on the dynamics of prices of the stocks considered.

In addition to this, also the average sentiment value of the news of the day
as well as the number of trades of the stock also influence the fluctuations
of prices in question. Moreover, it was also possible to observe that, among
the variables with a 7-day lag, the total return (gross dividends) is a
significant independent variable and thus also that the data concerning a
specific week has an impact on the dynamic of prices of the following period.
As it was stated earlier, these are the variables that were taken into
consideration in order to implement the neural network.

The following part of the analysis, in fact, involved building the neural
networks. To this regard, the aforementioned variables were taken into
consideration as independent variables, while the dependent variable was the
column indicating whether the specific company considered was present or
absent at a specific date in the STOXX Europe 600 Index.

The entire timespan was thus divided into two parts, the first part
corresponding to 80% of the total data and the second part corresponding to
20% of the total data; the former part was used to train the model while the
latter was used to test the model in question.

As it was stated earlier, the length of the timespan depended on the
availability of the data and therefore for some companies the model was built
considering approximately 450 periods, for example in the case of Cancom, and
for other companies the available periods were about 500 or more, as in the
case of Rotork PLC.

##  4 Outcome of data analysis

As it was stated earlier, the first part of the analysis involved the
employment of the logit model, through which it was possible to determine
which are the most significant independent variables with regard to the
analysis of the fluctuations of the prices; these variables were then selected
in order to build the neural networks for each company considered. From the
logit model it was thus possible to observe that among the most significant
variables there is the independent variable indicating whether the company at
the time instant considered is a constituent or is not part of the STOXX
Europe 600 Index. In other words, the dichotomic independent variable which
indicates with number 1 the fact that the company in question is part of the
Index and with number 0 the fact that the company considered is not part of
the STOXX Europe 600 Index is among the factors that determine the rise or the
decline in the prices of the specific stock. Similarly, it was possible to
observe that also the number of trades affects the dynamics of the
fluctuations of the prices of the stock. In addition to this, the third
independent variable that was identified as one of the variables that has an
impact on the prices of the stock is the average sentiment value of the news
of the day, which basically records the sentiment which can be derived from
the analysis of the information taken from the news. Finally, the only
variable belonging to the group of the independent variables with a 7-day lag
that was selected among the most significant variables is the total return
(gross dividends). More specifically, in this case the data concerning the
previous week was recorded and taken into consideration and therefore it was
thus possible to observe that also factors belonging to a different week
actually had an impact on the fluctuation of the following period of time.
These four factors, as it was mentioned earlier, were the four independent
variables that were taken into account to implement the neural networks for
the companies considered. Through the code written in Python it was possible
to use 80% of the data to train the model and the remaining 20% to test the
model in question. Therefore, this procedure was repeated for each company
taken into consideration and the various levels of accuracy were thus
determined. The levels of accuracy varied according to the specific company
considered, but overall the models built for the various companies showed
accuracy values of about 70%, as it can be seen from Table 1 (for the sake of
simplicity just 10 of the 50 companies are included in the following table).

Company Name | Accuracy
---|---
Cancom | 68,76%
Dechra Pharmaceuticals | 70,24%
Logitech International | 71,46%
Clariant | 69,49%
Deutsche Wohnen | 70,57%
DS Smith | 68,23%
Rotork | 69,65%
Roche Holding | 70,71%
RSA Insurance Group | 71,43%
RELX | 70,14%
Table 1: Values of accuracy

The authors, as it is specified also in the following section, believe that
there is still room for improvement, since the values of accuracy could
possibly increase if the model, for example, is constructed by taking into
consideration other combinations of significant independent variables.
Finally, the values of accuracy therefore provide a sufficiently satisfactory
result, since the values in question are above 65% and since they are
consistent throughout the various companies that have been taken into
consideration during the various parts of the analysis.

##  5 Conclusion

Through the thorough data analysis it was thus possible to observe and
understand the main determinants of the dynamics of the price. The first part
involved the careful and attentive data selection and cleaning processes.
Through this procedure it was possible to select with pre-set criteria which
companies had to be taken into consideration. In the following part of the
analysis the most significant independent variables were determined through
the employment of the logit model. These variables were then taken into
account to create the neural networks, which, applied to the different
companies considered, showed accuracy values of approximately the 70%, thus
providing a sufficiently satisfactory result.

Different types of future research could be related to this work. To this
regard, in fact, future work could deal with the expansion of the databases
considered, for example by taking into account more companies and longer
periods of time (when possible).

In addition to this, the focus of future work could also be quantifying the
effect that the inclusion or the exclusion from the STOXX Europe 600 Index has
on the fluctuations of prices.

Moreover, future work could also include the application of these analysis and
models to other indexes and therefore, in other words, to different economic
contexts.

In conclusion, a secondary type of analysis could also include the assessment
of the impact of these models on other aspects of these companies or
eventually indexes considered and therefore, on a larger scale, on various
economic and financial scenarios.

## References

  * [1] Ali, S. S., Mubeen, M., Lal, I., Hussain, A. (2018). Prediction of Stock Performance by Using Logistic Regression Model: Evidence from Pakistan Stock Exchange (PSX). _Asian Journal of Empirical Research_ , (7), 247?258.
  * [2] Belhouchet, S., Chouaibi, J. (2022). Do audit committee attributes improve integrated reporting quality? Evidence from European Companies listed on STOXX EUROPE 600 INDEX. _Research Square_ , pre print.
  * [3] Bonsón, E., Bednárová, M. (2015). CSR reporting practices of Eurozone companies. _Revista de Contabilidad_ 18 (2),182-193.
  * [4] Borucka, A. (2020). Logistic regression in modeling and assessment of transport services. _Open Engineering_ , 10(1), 26-34.
  * [5] Celine, S., Dominic, M. M, Devi, M. S. (2020). Logistic Regression for Employability Prediction. _International Journal of Innovative Technology and Exploring Engineering_ , 9(3).
  * [6] Chen, Y., Mantegna, R. N., Pantelous, A. A., Zuev, K. M. (2018). A dynamic analysis of S&P 500, FTSE 100 and EURO STOXX 50 indices under different exchange rates. _PLoS ONE_ 13(3): e0194067.
  * [7] Christiansen, H., Koldertsova, A. (2009). The role of stock exchanges in corporate governance. _Financial Market Trends_ , 1-32.
  * [8] De Blasis, R., Petroni, F. (2021). Price Leadership and Volatility Linkages between Oil and Renewable Energy Firms during the COVID-19 Pandemic. _Energies_ , 14, 2608.
  * [9] Giannarakis, G., Partalidou, X., Zafeiriou, E., Sariannidis, N. (2016). An analysis of United States on Dow Jones Sustainability Index. _Investment Management and Financial Innovations_ , 13(3-2), 353-361.
  * [10] Kyriazis, A. N. (2021). Investigating the nexus between European major and sectoral stock indices, gold and oil during the COVID 19 pandemic. _SN Business & Economic_, 1, 57.
  * [11] Lachiheb, O., Gouider, M. S. (2018). A hierarchical Deep neural network design for stock returns prediction. _Procedia Computer Science_ 126, 264-272.
  * [12] Langenstein, L., Uzík, M., Warias, R. (2021). Winner Strategies in Crisis. SHS Web of Conferences 92. Globalization and its Socio-Economic Consequences 2020.
  * [13] Lee, Y. H., Hao, F. (2012). Oil and S&P 500 Markets: Evidence from the Nonlinear Model. _International Journal of Economics and Financial Issues_ , Econjournals, 2(3), 272-280.
  * [14] Lento, C., Gradojevic, N. (2021). S&P 500 Index Price Spillovers around the COVID-19 Market Meltdown. _Journal of Risk and Financial Management_ 14, 330.
  * [15] Mehtab, S., Sen, J (2019). A Robust Predictive Model for Stock Price Prediction Using Deep Learning and Natural Language Processing. Available at SSRN:https://ssrn.com/abstract=3502624 or http://dx.doi.org/10.2139/ssrn.3502624.
  * [16] Metghalchi, M., Chen, C. P., Hajilee, M. (2016). Moving average trading rules for NASDAQ composite index. _Applied Finance Letters_ , 5(2), 45-57.
  * [17] Mironiuc, M., Robu, M. A. (2013). Obtaining a Practical Model for Estimating Stock Performance on an Emerging Market Using Logistic Regression Analysis. _Procedia - Social and Behavioral Sciences_ 81, 422-427.
  * [18] Moghaddam, A. H., Moghaddam, M. H., Esfandyari, M. (2016). Stock market index prediction using artificial neural network. _Journal of Economics, Finance and Administrative Science_ , 21(41), 89-93.
  * [19] Salvi, A., Teti, E., Giakoumelou, A. (2021). The Relationship between LBOs, IPOs and Market Risk Premium: An Empirical Analysis of the European Market. _International Journal of Business and Management_ , Canadian Center of Science and Education, 14(7), 160-160.
  * [20] Shahvaroughi Farahani, M., Razavi Hajiagha, S. H. (2021). Forecasting stock price using integrated artificial neural network and metaheuristic algorithms compared to time series models. _Soft Comput_ , 25, 8483-8513.
  * [21] Smita, M. (2021). Logistic Regression Model For Predicting Performance of S&P BSE30 Company Using IBM SPSS. _International Journal of Mathematics Trends and Technology_ , 67(7), 118-134.
  * [22] Song, Y. G., Zhou, Y. L., Han, R. J. (2018). Neural networks for stock price prediction. _Journal of Difference Equations and Applications_ , 1?14.
  * [23] Zaini, B. J., Mansor, R., Yusof, N., Sang, B. H. (2019). Classify Stock Market Movement Based on Technical Analysis Indicators Using Logistic Regression. _Journal of Advanced Research in Business and Management Studies_ , 14 (1), 35-41.


# Semiparametric Inference for Regression-Discontinuity Designs

Rong J.B. Zhu and Weiwei Jiang
Institute of Science and Technology for Brain-inspired Intelligence
Fudan University, Shangha China
E-mail<EMAIL_ADDRESS>

###### Abstract

Treatment effects in regression discontinuity designs (RDDs) are often
estimated using local regression methods. However, global approximation
methods are generally deemed inefficient. In this paper, we propose a
semiparametric framework tailored for estimating treatment effects in RDDs.
Our global approach conceptualizes the identification of treatment effects
within RDDs as a partially linear modeling problem, with the linear component
capturing the treatment effect. Furthermore, we utilize the $P$-spline method
to approximate the nonparametric function and develop procedures for inferring
treatment effects within this framework. We demonstrate through Monte Carlo
simulations that the proposed method performs well across various scenarios.
Furthermore, we illustrate using real-world datasets that our global approach
may result in more reliable inference.

Keywords: Causal inference; Regression discontinuity design; Partially linear
model; $P$-spline.
JEL Classification: C13, C21, C23.

##  1 Introduction

The regression discontinuity design (RDD) was initially introduced by
Thistlewaite and Campbell (1960) to examine the influence of merit awards on
students’ the future academic outcomes. Since then, it has evolved into one of
the most reputable strategies for estimating treatment effects across
disciplines, such as economics, political science, and medicine. In an RDD,
units receive scores, and treatment allocation is determined by whether the
score exceeds a predetermined cutoff value: units scoring above the cutoff are
designated to the treatment condition, while those scoring below are assigned
to the control condition. This treatment assignment rule creates a
discontinuity in the probability of receiving treatment, allowing researchers
to learn about the treatment effect by comparing units barely above and barely
below the cutoff. For comprehensive reviews of the RDD literature, refer to
Imbens and Lemieux (2008) and Cattaneo et al. (2019).

The local linear or quadratic approximation (Fan, 1992; Fan and Gijbels, 1996)
is a commonly employed approach in RDDs (Hahn et al., 2001). Specifically,
researchers typically estimate the treatment effect locally around the cutoff
point, which depends on the chosen bandwidth. Studies have investigated
optimal bandwidth selection for local linear and quadratic regression
specifications (Imbens and Kalyanaraman, 2012; Calonico et al., 2014, 2015).
Another approach that researchers have explored is the utilization of a global
polynomial regression approach (Lee and Lemieux, 2010). In practice,
researchers often employ high order (up to fifth or sixth-order) polynomials,
choosing the degree of the polynomial using statistical information criteria
or cross-validation. However, Gelman and Imbens (2019) have argued against the
use of high-order global polynomial approximations, advocating instead for
inference based on local low-order polynomials, such as local linear or local
quadratic models.

In this article, we delve into a semiparametric inference framework for RDDs,
proposing a global estimation method, which could serve as an alternative to
local estimation methods. Our primary focus is on identifying the treatment
effect within the framework of a partially linear model. This entails
capturing the treatment effect through the coefficient of the treatment
assignment variable in the linear component of the partially linear model. By
employing $P$-spline to approximate the nonparameteric part, we reformulate
the semiparametric RDD framework into a linear mixed model regression.
Subsequently, we derive the estimate of the treatment effect within the linear
mixed model regression framework. Theoretically, we establish the asymptotic
normality of the estimate under regular conditions, providing an inference
procedure based on this normality. Empirically, our experiments on simulated
datasets support the assertion that our global estimation provides a viable
alternative to local estimation methods. Furthermore, our examination of real-
world datasets demonstrates that our global estimation approach may yield more
reliable inferences.

The rest of the article proceeds as follows. Section 2 proposes a
semiparametric framework that incorporates the partially linear model to
estimate treatment effects. Section 3 discusses estimation and inference
procedures. Section 4 presents the empirical analysis. Section 5 showcases the
performance on real datasets. Concluding remarks are provided in Section 6.
Proofs are dispersed between the main text and the appendix.

##  2 A Semiparametric framework for regression-discontinuity designs

###  2.1 Setup and assumptions

Employing the potential outcome approach, let $(Y_{i}(0),Y_{i}(1))$ denote the
pair of potential outcomes for individual $i=1,\cdots,n$ with
$w_{i}\in\\{0,1\\}$ representing the treatment received. There are $n_{1}$
units that receive the treatment, while the remaining $n_{0}=n-n_{1}$ units
serve as controls. Each individual $i$ is assigned only one treatment, denoted
by $w_{i}$, resulting in the observation of a single realized outcome,
$y_{i}=Y_{i}(w_{i})$.

We study the RDDs, particularly the sharp case, where the treatment is a
monotone, deterministic function of a pretreatment variable (forcing variable)
$x_{i}$ being on either side of a fixed and known threshold $x_{0}$, i.e.,
$w_{i}=\mathbb{I}\\{x_{i}\geq x_{0}\\}$. Denote the conditional expectations
of the potential outcomes as

$\displaystyle\mu_{1}(x)=$ $\displaystyle\mathbb{E}[Y_{i}(1)|x_{i}=x],$
$\displaystyle\mu_{0}(x)=$ $\displaystyle\mathbb{E}[Y_{i}(0)|x_{i}=x].$

The average effect at the threshold $x_{0}$ is defined as

$\tau_{x_{0}}=\mu_{1}(x_{0})-\mu_{0}(x_{0}).$

We impose the following assumptions:

######  Assumption 1.

Assume that the potential outcomes are modeled as

$\displaystyle Y_{i}(1)=$ $\displaystyle\mu_{1}(x_{i})+\epsilon_{i}(1),$
$\displaystyle Y_{i}(0)=$ $\displaystyle\mu_{0}(x_{i})+\epsilon_{i}(0).$

######  Assumption 2.

Assume that the mean functions $\mu_{1}(x)$ and $\mu_{0}(x)$ are continuous.

######  Assumption 3.

Assume that the model errors $(\epsilon_{i}(1),\epsilon_{i}(0))_{i=1}^{n}$
satisfy the following properties:
$\\{(\epsilon_{1}(0),\epsilon_{1}(1))^{\top},\cdots,(\epsilon_{n}(0),\epsilon_{n}(1))^{\top}\\}$
is an independent and identically distributed sequence with
$(\epsilon_{i}(0),\epsilon_{i}(1))^{\top}\sim\mathcal{N}(0,\sigma^{2}\mathbf{I}_{2})$.

Assumptions 1 & 2 are standard in RDDs (Imbens and Lemieux, 2008). In
Assumption 3, for convenience, we assume homoscedasticity in $\epsilon_{i}(1)$
and $\epsilon_{i}(0)$. Nevertheless, this assumption can be extended to
accommodate heteroscedasticity in these errors, thereby allowing for two
different variances.

###  2.2 Semiparametric RDD framework

Now, we introduce a semiparametric framework for estimating the average effect
of the treatment. Denote
$\epsilon_{i}=\epsilon_{i}(1)w_{i}+\epsilon_{i}(0)(1-w_{i})$, and define

$\displaystyle
g(x_{i})=\mu_{0}(x_{i})\mathbb{I}\\{x_{i}<x_{0}\\}+(\mu_{1}(x_{i})-\tau_{x_{0}})\mathbb{I}\\{x_{i}\geq
x_{0}\\}.$ (1)

######  Theorem 1.

Suppose that Assumptions 1, 2, and 3 hold. We have that the observations are
generated from the following model:

$\displaystyle y_{i}=$ $\displaystyle\tau_{x_{0}}w_{i}+g(x_{i})+\epsilon_{i}.$
(2)

Moreover, $g(x)$ is continuous at $x_{0}$, and
$\epsilon_{1},\cdots,\epsilon_{n}$ are independent with
$\text{Var}(\epsilon_{i})=\sigma^{2}$.

###### Proof.

Model (2) is a direct result from Assumption 1 and the definition of
$g(x_{i})$ in (1).

From the definition of $g(x)$,

$\displaystyle g(x_{0})=$
$\displaystyle\mu_{1}(x_{0})-\tau_{x_{0}}=\mu_{0}(x_{0})$
$\displaystyle\lim_{x\rightarrow x_{0}^{+}}g(x)=$
$\displaystyle\lim_{x\rightarrow
x_{0}^{+}}\mu_{1}(x)-\tau_{x_{0}}=\mu_{0}(x_{0})$
$\displaystyle\lim_{x\rightarrow x_{0}^{-}}g(x)=$
$\displaystyle\lim_{x\rightarrow x_{0}^{+}}\mu_{0}(x)=\mu_{0}(x_{0}),$

where the last steps in the last two equations above are from Assumption 2.
The three equations above imply that $g(x)$ is continuous at $x_{0}$.

From Assumption 3, we have that $\epsilon_{1},\cdots,\epsilon_{n}$ are
independent, and

$\displaystyle\text{Var}(\epsilon_{i})=$
$\displaystyle\text{Var}[\epsilon_{i}(1)w_{i}+\epsilon_{i}(0)(1-w_{i})]=\sigma^{2}.$

Therefore, Theorem 1 is proved. ∎

Theorem 1 establishes the identifiability of the treatment effect for RDDs
using the partially linear model (2), wherein inference on the treatment
effect can be construed as inference on the coefficient of $w_{i}$.

##  3 Semiparametric inference

###  3.1 $P$-spline approach on modeling $g(x)$

We approximate $g(x)$ using $P$-splines with the following function:

$\displaystyle g_{m}(x)=$
$\displaystyle\beta_{0}+\beta_{1}x+\cdots\beta_{p}x^{p}+\sum\limits_{k=1}^{K}\gamma_{k}(x-\kappa_{k})_{+}^{p},$
(3)

where $p$ (taking values of 1 or 2) represents the degree of the spline,
$(x)_{+}^{p}=x^{p}\mathbb{I}\\{x\geq 0\\}$, and $\kappa_{1}<\cdots<\kappa_{K}$
is a set of knots over the range of $x$. Given the knots
$\boldsymbol{\kappa}=(\kappa_{1},\cdots,\kappa_{K})^{\top}$ and $p$, Model (3)
use the truncated polynomial spline basis
$\\{1,x,\cdots,x^{p},(x-\kappa_{1})_{+}^{p},\cdots,(x-\kappa_{K})_{+}^{p}\\}$
to approximate the function $g(x)$. Denote
$\mathbf{x}=(1,x,\cdots,x^{p})^{\top}$,
$\mathbf{z}=((x-\kappa_{1})_{+}^{p},\cdots,(x-\kappa_{K})_{+}^{p})^{\top}$,
$\boldsymbol{\beta}=(\beta_{0},\cdots,\beta_{p})^{\top}$, and
$\boldsymbol{\gamma}=(\gamma_{1},\cdots,\gamma_{K})^{\top}$. Model (3) is
rewritten as

$\displaystyle g_{m}(x)=$
$\displaystyle\mathbf{x}^{\top}\boldsymbol{\beta}+\mathbf{z}^{\top}\boldsymbol{\gamma}.$
(4)

Provided that the knots are sufficiently spread out and the number of knots
$K$ is sufficiently large, the approximation above is of a high level of
accuracy. The model may encounter overfitting issues due to potential over-
parameterization. To mitigate this issue, a penalty is imposed on the
magnitude of the spline parameters $\boldsymbol{\gamma}$, treating
$\gamma_{k}\sim N(0,\sigma_{\gamma}^{2})$ as a random effect with variance
$\sigma_{\gamma}^{2}$ (Ruppert et al., 2003). Empirical studies by Ruppert
(2002) suggested that the approximation bias using the $P$-spline approach is
negligible in univariate nonparametric regression cases. Theoretical analysis
by Claeskens et al. (2009) demonstrated that the bias is
$O(n^{-(p+1)/(2p+1)})$ assuming $\lambda=O(1)$ and $K=O(n^{1/(2p+3)})$, which
is negligible.

Ruppert et al. (2003) provided recommendations on knot selection. The knots
are frequently placed at equally spaced quantiles of the distribution of $x$,
with $K$ chosen to be large relative to the size of the dataset. For
univariate $x$, a typical knot choice is to have one knot every four or five
observations, with a maximum number ranging from 35 to 50.

The family $\\{g_{m}(x)\\}$ includes functions with discontinuous derivatives
or discontinuities themselves. Such occurrences naturally arise in splines
when a knot is placed at the discontinuity location. Consequently, our
proposed method can effectively handle functions $g$ with isolated
discontinuities or discontinuous derivatives.

###  3.2 Reformulation to linear mixed model

As the $P$-spline approach can be viewed as random-effects models (Ruppert,
2002), it is natural to reformulate the semiparametric RDD framework to linear
mixed model regression. Define $\mathbf{y}=(y_{1},\cdots,y_{n})^{\top}$,
$\mathbf{w}=(w_{1},\cdots,w_{n})^{\top}$,
$\mathbf{X}=(\mathbf{x}_{1},\cdots,\mathbf{x}_{n})^{\top}$,
$\mathbf{Z}=(\mathbf{z}_{1},\cdots,\mathbf{z}_{n})^{\top}$, and
$\boldsymbol{\epsilon}=(\epsilon_{1},\cdots,\epsilon_{n})^{\top}$. We treat
the data as if they came from the regression model

$\displaystyle\mathbf{y}=\tau_{x_{0}}\mathbf{w}+\mathbf{X}\boldsymbol{\beta}+\mathbf{Z}\boldsymbol{\gamma}+\boldsymbol{\epsilon},$
(5)

where
$\boldsymbol{\gamma}\sim\mathcal{N}(0,\sigma_{\gamma}^{2}\mathbf{I}_{K})$ and
$\boldsymbol{\epsilon}\sim\mathcal{N}(0,\sigma^{2}\mathbf{I}_{n})$.

Denoting
$\boldsymbol{\theta}=(\tau_{x_{0}},\boldsymbol{\beta}^{\top})^{\top}$, we
estimate the coefficients $\boldsymbol{\theta}$ using generalized least
squares. In the model, the covariance matrix of $\mathbf{y}$ is given by

$\mathbf{V}=\sigma_{\gamma}^{2}\mathbf{Z}\mathbf{Z}^{\top}+\sigma^{2}\mathbf{I}_{n}.$

We have $\mathbf{V}^{-1}=\sigma^{-2}\mathbf{Q}_{z}$, where

$\mathbf{Q}_{z}=\mathbf{I}-\mathbf{Z}(\mathbf{Z}^{\top}\mathbf{Z}+\sigma^{2}\sigma_{\gamma}^{-2}\mathbf{I}_{K})^{-1}\mathbf{Z}^{\top}.$

Denote $\mathbf{U}=(\mathbf{w},\mathbf{X})$. Assuming known variances
$\sigma_{\gamma}^{2}$ and $\sigma^{2}$, we estimate $\boldsymbol{\theta}$
using the generalized least squares estimate:

$\displaystyle\tilde{\boldsymbol{\theta}}=$
$\displaystyle(\mathbf{U}^{\top}\mathbf{V}^{-1}\mathbf{U})^{-1}\mathbf{U}^{\top}\mathbf{V}^{-1}\mathbf{y}=\left(\mathbf{U}^{\top}\mathbf{Q}_{z}\mathbf{U}\right)^{-1}\left(\mathbf{U}^{\top}\mathbf{Q}_{z}\mathbf{y}\right)$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\mathbf{w}^{\top}\mathbf{Q}_{z}\mathbf{w}&\mathbf{w}^{\top}\mathbf{Q}_{z}\mathbf{X}\\\
\mathbf{X}^{\top}\mathbf{Q}_{z}\mathbf{w}&\mathbf{X}^{\top}\mathbf{Q}_{z}\mathbf{X}\end{array}\right)^{-1}\left(\begin{array}[]{c}\mathbf{w}^{\top}\mathbf{Q}_{z}\mathbf{y}\\\
\mathbf{X}^{\top}\mathbf{Q}_{z}\mathbf{y}\end{array}\right).$ (10)

The estimates of $\sigma_{\gamma}^{2}$ and $\sigma^{2}$ can be obtained using
maximum likelihood or restricted maximum likelihood (Patterson and Thompson,
1971), which are implemented by the “nlme” (Pinheiro et al., 2023) or “lme4”
(Bates et al., 2015) packages in R. Once consistent estimators
$\hat{\sigma}_{\gamma}^{2}$ and $\hat{\sigma}^{2}$ for $\sigma_{\gamma}^{2}$
and $\sigma^{2}$ are obtained, they can be inserted into (3.2) to obtain the
feasible GLS estimator.

$\displaystyle\hat{\boldsymbol{\theta}}=$
$\displaystyle\left(\begin{array}[]{cc}\mathbf{w}^{\top}\hat{\mathbf{Q}}_{z}\mathbf{w}&\mathbf{w}^{\top}\hat{\mathbf{Q}}_{z}\mathbf{X}\\\
\mathbf{X}^{\top}\hat{\mathbf{Q}}_{z}\mathbf{w}&\mathbf{X}^{\top}\hat{\mathbf{Q}}_{z}\mathbf{X}\end{array}\right)^{-1}\left(\begin{array}[]{c}\mathbf{w}^{\top}\hat{\mathbf{Q}}_{z}\mathbf{y}\\\
\mathbf{X}^{\top}\hat{\mathbf{Q}}_{z}\mathbf{y}\end{array}\right),$ (15)

where
$\hat{\mathbf{Q}}_{z}=\mathbf{I}-\mathbf{Z}(\mathbf{Z}^{\top}\mathbf{Z}+\hat{\sigma}^{2}\hat{\sigma}_{\gamma}^{-2}\mathbf{I}_{K})^{-1}\mathbf{Z}^{\top}$.

###  3.3 Inference

Utilizing standard results on asymptotic normality of maximum likelihood
estimates (MLE) in mixed effects models, we present the asymptotic normality
of $\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}$ in the following theorem.

Denote $\rho_{1}=\lim_{n\rightarrow\infty}n_{1}/n$,
$\boldsymbol{\mu}_{z}=\lim_{n\rightarrow\infty}n^{-1}\sum\limits_{i=1}^{n}\mathbf{z}_{i}$,
and
$\boldsymbol{\Sigma}_{z}=\lim_{n\rightarrow\infty}n^{-1}\sum\limits_{i=1}^{n}\mathbf{z}_{i}\mathbf{z}_{i}^{\top}$.
We impose the following assumption for the theorem.

######  Assumption 4.

(a) $K\rightarrow\infty$ as $n\rightarrow\infty$.
(b) $\lim_{n\rightarrow\infty}K/n=0$.
(c) $r=\lim_{n\rightarrow\infty}\text{rank}(\mathbf{Z})/K>0$.
(d) $\rho_{1}$ satisfies the following inequalities and $\rho_{1}>0$ and
$\rho_{1}<1/(1+\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z})$.

Parts (a) & (b) in Assumption 4 impose constraints on the number of knots $K$.
These conditions require that the number of knots approaches infinity as
$n\rightarrow\infty$, and that the number of knots is of a smaller order than
$n$ (i.e., $o(n)$). Part (c) permits the rank of the design matrix of
$\mathbf{Z}$ to be smaller than $K$ but of the same order.

Part (d) imposes the conditions on the ratio of treated units. Specifically,
$\rho_{1}>0$ requires the presence of treated units in the vicinity of the
threshold. Now, our attention shifts to the condition
$\rho_{1}<1/(1+\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z})$.
Denoting $\mathbf{V}_{z}$ as the variance of $\mathbf{z}_{i}$, we have that

$\displaystyle\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}=$
$\displaystyle\boldsymbol{\mu}_{z}^{\top}(\boldsymbol{\mu}_{z}\boldsymbol{\mu}_{z}^{\top}+\mathbf{V}_{z})^{-1}\boldsymbol{\mu}_{z}=\frac{\boldsymbol{\mu}_{z}^{\top}\mathbf{V}_{z}^{-1}\boldsymbol{\mu}_{z}}{1+\boldsymbol{\mu}_{z}^{\top}\mathbf{V}_{z}^{-1}\boldsymbol{\mu}_{z}}.$

This leads to the expression for the condition:

$\rho_{1}<(1+\boldsymbol{\mu}_{z}^{\top}\mathbf{V}_{z}^{-1}\boldsymbol{\mu}_{z})/(1+2\boldsymbol{\mu}_{z}^{\top}\mathbf{V}_{z}^{-1}\boldsymbol{\mu}_{z}).$

As
$\boldsymbol{\mu}_{z}^{\top}\mathbf{V}_{z}^{-1}\boldsymbol{\mu}_{z}\rightarrow
0$, the requirement simplifies to $\rho_{1}<1$. Conversely, as
$\boldsymbol{\mu}_{z}^{\top}\mathbf{V}_{z}^{-1}\boldsymbol{\mu}_{z}\rightarrow\infty$,
the requirement becomes $\rho_{1}<1/2$. This indicates that as
$\boldsymbol{\mu}_{z}^{\top}\mathbf{V}_{z}^{-1}\boldsymbol{\mu}_{z}$
decreases, the requirement on the ratio of treated units becomes less strict.
The term $\boldsymbol{\mu}_{z}^{\top}\mathbf{V}_{z}^{-1}\boldsymbol{\mu}_{z}$
quantifies the dispersion of the data $\mathbf{z}_{i}$, $i=1,\cdots,n$, within
the space defined by the inverse of the variance-covariance matrix
$\mathbf{V}_{z}^{-1}$. A larger value of
$\boldsymbol{\mu}_{z}^{\top}\mathbf{V}_{z}^{-1}\boldsymbol{\mu}_{z}$ indicates
that the data points tend to cluster closer to the mean vector, suggesting
less dispersion or variability in the data.

######  Theorem 2.

Define $(p+2)\times(p+2)$ matrix $\mathbf{J}$ given by

$\mathbf{J}=\lim_{n\rightarrow\infty}\mathbf{U}^{\top}\mathbf{V}^{-1}\mathbf{U}/n.$

Under the model (5), assuming that Assumption 4 holds, we have that

$\displaystyle\sqrt{n}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta})\rightarrow_{d}\mathcal{N}(0,\mathbf{J}^{-1}),\text{
as }n\rightarrow\infty.$ (16)

###### Proof.

The result in (16) from an application of Miller (1977). See details in
Appendix A. ∎

Theorem 2 is proven in the case where the unknown variance components are
substituted by the MLE estimates. However, in the scenario where the unknown
variance components are estimated using restricted maximum likelihood
estimation, asymptotic normality still holds, as demonstrated by Jiang (1998).
Nonetheless, further elaboration is required for a complete derivation,
applying the result outlined in Jiang (1998).

The primary outcome of Theorem 2 is to offer inference on $\tau_{x_{0}}$.
Noting
$\hat{\boldsymbol{\theta}}=(\hat{\tau}_{x_{0}},\hat{\boldsymbol{\beta}}^{\top})^{\top}$,
the asymptotic variance of $\hat{\tau}_{x_{0}}$ is obtained from
$(\mathbf{U}^{\top}\mathbf{V}^{-1}\mathbf{U})^{-1}$ using the block-wise
inversion formula. Applying the block-wise inversion formula to
$(\mathbf{U}^{\top}\mathbf{V}^{-1}\mathbf{U})^{-1}$ implies that its
$(1,1)$-th element is given by

$V_{\tau}=\left(\mathbf{w}^{\top}\mathbf{Q}_{z}\mathbf{w}-\mathbf{w}^{\top}\mathbf{Q}_{z}\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}_{z}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Q}_{z}\mathbf{w}\right)^{-1}.$

Therefore, Theorem 2 implies:

$\displaystyle
V_{\tau}^{-1/2}(\hat{\tau}_{x_{0}}-\tau_{x_{0}})\rightarrow_{d}\mathcal{N}(0,1),\text{
as }n\rightarrow\infty.$ (17)

The variance $V_{\tau}$ is estimated as

$\hat{V}_{\tau}=\left(\mathbf{w}^{\top}\hat{\mathbf{Q}}_{z}\mathbf{w}-\mathbf{w}^{\top}\hat{\mathbf{Q}}_{z}\mathbf{X}(\mathbf{X}^{\top}\hat{\mathbf{Q}}_{z}\mathbf{X})^{-1}\mathbf{X}^{\top}\hat{\mathbf{Q}}_{z}\mathbf{w}\right)^{-1}.$

Hence, an asymptotic point-wise $100(1-\alpha)\%$ confidence interval for
$\tau_{x_{0}}$ is given by

$\displaystyle\left(\hat{\tau}_{x_{0}}-z_{\alpha/2}\hat{V}_{\tau}^{1/2},\hat{\tau}_{x_{0}}+z_{\alpha/2}\hat{V}_{\tau}^{1/2}\right),$
(18)

where $z_{\alpha/2}$ is the upper $\alpha/2$-quantile from a standard Gaussian
distribution.

##  4 Numerical experiments

In this section, we conduct Monte Carlo simulations to examine the finite
sample performance of our proposed method. We carry out 5000 replications,
generating a random sample $\\{(x_{i},\epsilon_{i})^{\top}:i=1,\cdots,n\\}$
with a sample size of $n=50$ or 200 for each replication. The covariate
$x_{i}=2z_{i}-1$, where $z_{i}$ follows a beta distribution of
$\text{beta}(3,3)$, and the error term $\epsilon_{i}$ follows a normal
distribution $\mathcal{N}(0,\sigma^{2})$ with $\sigma=0.2$.

In our method’s implementation, we set the number of knots $K=20$. We also
explore the automatic way of choosing $K$ as proposed by Ruppert et al.
(2003), yielding similar results. Once $K$ is determined, the knots are placed
at fixed quantiles of $\\{x_{1},\cdots,x_{n}\\}$. We consider two values for
$p=1,2$, denoted as “PL1” and “PL2”, respectively.

We consider three regression functions labeled M1, M2, and M3, respectively:

M1: $\displaystyle\mu(x)=x\text{ for }x<0;\ \ \mu(x)=\tau+x\text{ for }x\geq
0$ M2: $\displaystyle\mu(x)=1-\exp(-x)\text{ for }x<0;\ \
\mu(x)=\exp(x)-1\text{ for }x\geq 0$ M3:
$\displaystyle\mu(x)=0.48+1.27x+7.18x^{2}+20.21x^{3}+21.54x^{4}+7.33x^{5}\text{
for }x<0;$
$\displaystyle\mu(x)=\tau+0.48+0.84x-3x^{2}+7.99x^{3}-9.01x^{4}+3.56x^{5}\text{
for }x\geq 0.$

In M1, there exists a consistent linear trend between $\mu_{0}(x)$ and
$\mu_{1}(x)$. In M2, both $\mu_{0}(x)$ and $\mu_{1}(x)$ exhibit nonlinearity.
M3 is from Lee (2008). The outcomes are generated as
$y_{i}=\mu(x_{i})+\epsilon_{i}$, $i=1,\cdots,n$, for each model. We compare
our approach to the local linear RD estimator ($p=1$), denoted as “IK1”
(Imbens and Kalyanaraman, 2012). Additionally, we implemented local quadratic
RD estimator ($p=2$), but observed that it performed worse than the local
linear estimator. Consequently, we opted not to report the performance of the
local quadratic RD estimator.

Table 1: The performance of the PL method under various model settings.
“RMSE”: root mean square error, “Bias”: bias values, “EC”: empirical coverage
rate, and “IL”: average confidence interval Length.

$n$ |  | M1 | M2 | M3
---|---|---|---|---
RMSE | Bias | EC | IL | RMSE | Bias | EC | IL | RMSE | Bias | EC | IL
50 | IK1 | 0.127 | 0.005 | 0.947 | 0.465 | 0.146 | -0.043 | 0.894 | 0.488 | 0.135 | 0.037 | 0.926 | 0.472
PL1 | 0.111 | 0.003 | 0.945 | 0.425 | 0.155 | -0.087 | 0.833 | 0.448 | 0.116 | 0.033 | 0.943 | 0.427
PL2 | 0.120 | 0.005 | 0.942 | 0.447 | 0.144 | -0.050 | 0.882 | 0.472 | 0.126 | 0.036 | 0.937 | 0.449
200 | IK1 | 0.059 | -0.001 | 0.937 | 0.216 | 0.083 | -0.036 | 0.848 | 0.251 | 0.066 | 0.027 | 0.913 | 0.218
PL1 | 0.053 | -0.001 | 0.950 | 0.204 | 0.102 | -0.060 | 0.700 | 0.254 | 0.061 | 0.029 | 0.908 | 0.207
PL2 | 0.057 | -0.001 | 0.945 | 0.211 | 0.080 | -0.033 | 0.870 | 0.254 | 0.065 | 0.029 | 0.910 | 0.213

Table 1 presents simulation results under $\tau=0$, featuring metrics such as
root-mean-square error (“RMSE”), bias (“Bias”), empirical coverage rate of the
95% CI (“EC”), and the average empirical interval length (“IL”). The results
demonstrate a reduction of 10% in RMSE for M1 and M3 with PL1 compared to IK1,
alongside no noticeable increases in bias. Moreover, there are shorter ILs,
with EC approaching 95%. Empirically, the bias resulting from the $P$-spline
approximation is found to be negligible, as discussed in Section 3.1. In the
case of M2, PL2 outperforms IK1, while PL1 performs worse. This suggests that
employing $p=2$ for nonlinear functions may enhance the effectiveness of the
PL method. Additionally, a comparison of the performance under $n=50$ and
$n=200$ reveals consistent performance across sample sizes, indicating its
suitability for various applications.

Figure 1: An example model (M2) in the simulation setting: the fitted
functions.

The performance on inference. Apart from the estimation performance, we also
investigate the inference performance. Table 2 presents the type-i@ and type-
ii@ errors for each method. According to the table, PL1 exhibits a type-i@
error of 0.05, surpassing IK1 in M1 and showing a slightly higher type-i@
error than IK1 in M3. In the case of M2, PL2 outperforms IK1, while PL1
performs worse. For $\tau=0.1$, PL1 and PL2 have lower type-ii@ error than IK1
in M1 and M3, with PL2 also showcasing lower type-ii@ error than IK1 in M2. As
$\tau$ increases to 0.3, the PL method displays slightly higher type-ii@
error, though the difference is marginal. In conclusion, our approach offers a
viable alternative to the local estimator, delivering comparable performance
in inference,

Table 2: Performance on Testing

| M1 | M2 | M3 | M1 | M2 | M3 | M1 | M2 | M3
---|---|---|---|---|---|---|---|---|---
| Type-i@ error($\tau=0$) | Type-ii@ error($\tau=0.1$) | Type-ii@ error($\tau=0.3$)
IK1 | 0.063 | 0.152 | 0.087 | 0.545 | 0.820 | 0.379 | 0.001 | 0.015 | 0.000
PL1 | 0.050 | 0.300 | 0.092 | 0.536 | 0.872 | 0.316 | 0.003 | 0.038 | 0.002
PL2 | 0.055 | 0.130 | 0.090 | 0.537 | 0.810 | 0.344 | 0.002 | 0.016 | 0.001

Impact of the distribution of $x_{i}$. We investigate the performance of the
PL method under varying distributions of $x_{i}$. We generate $x_{i}=2z_{i}-1$
with $z_{i}\sim\text{Beta}(\alpha,\beta)$, considering two cases:
$(\alpha,\beta)=(0.5,0.5)$ and $(\alpha,\beta)=(4,4)$. It is notable that
there are more data points around the cut-off point in the case of
$(\alpha,\beta)=(4,4)$ compared to $(\alpha,\beta)=(0.5,0.5)$. Table 3
demonstrates the robustness of the PL method across different distributions.

Table 3: The performance of the PL method under various distributions of
$x_{i}$.

$(\alpha,\beta)$ |  | M1 | M2 | M3
---|---|---|---|---
RMSE | Bias | EC | IL | RMSE | Bias | EC | IL | RMSE | Bias | EC | IL
$(0.5,0.5)$ | IK1 | 0.077 | 0.003 | 0.925 | 0.272 | 0.134 | -0.069 | 0.832 | 0.375 | 0.092 | 0.017 | 0.935 | 0.356
PL1 | 0.070 | 0.003 | 0.944 | 0.265 | 0.136 | -0.071 | 0.868 | 0.417 | 0.110 | 0.039 | 0.943 | 0.438
PL2 | 0.079 | 0.002 | 0.929 | 0.281 | 0.112 | -0.044 | 0.912 | 0.388 | 0.109 | 0.011 | 0.943 | 0.418
$(4,4)$ | IK1 | 0.058 | -0.003 | 0.943 | 0.213 | 0.075 | -0.034 | 0.870 | 0.239 | 0.068 | 0.033 | 0.898 | 0.215
PL1 | 0.052 | -0.003 | 0.947 | 0.202 | 0.095 | -0.061 | 0.689 | 0.231 | 0.059 | 0.027 | 0.918 | 0.204
PL2 | 0.055 | -0.003 | 0.945 | 0.207 | 0.076 | -0.035 | 0.847 | 0.239 | 0.064 | 0.031 | 0.904 | 0.210

Impact of the noise. Next, we explore different values of $\sigma$,
specifically $\sigma=0.1,0.8$, and report the results in Table 4. Our method
exhibits robust performance under both settings, indicating its resilience
across different levels of noise variance.

Table 4: The performance of the PL method under various $\sigma$ values.

$\sigma$ |  | M1 | M2 | M3
---|---|---|---|---
RMSE | Bias | EC | IL | RMSE | Bias | EC | IL | RMSE | Bias | EC | IL
0.1 | IK1 | 0.029 | -0.001 | 0.939 | 0.107 | 0.048 | -0.026 | 0.843 | 0.137 | 0.042 | 0.028 | 0.808 | 0.110
PL1 | 0.027 | -0.001 | 0.950 | 0.102 | 0.046 | -0.019 | 0.915 | 0.166 | 0.043 | 0.031 | 0.787 | 0.109
PL2 | 0.029 | -0.001 | 0.945 | 0.106 | 0.040 | -0.017 | 0.911 | 0.140 | 0.043 | 0.030 | 0.796 | 0.111
0.8 | IK1 | 0.238 | -0.007 | 0.936 | 0.869 | 0.253 | -0.058 | 0.919 | 0.884 | 0.239 | 0.021 | 0.940 | 0.870
PL1 | 0.213 | -0.006 | 0.950 | 0.817 | 0.249 | -0.109 | 0.906 | 0.828 | 0.214 | 0.024 | 0.947 | 0.818
PL2 | 0.230 | -0.005 | 0.945 | 0.846 | 0.250 | -0.071 | 0.919 | 0.864 | 0.231 | 0.025 | 0.944 | 0.846

##  5 Real data

We illustrate our method here using two applications of RDDs. The first
dateset is sourced from U.S. House elections (Lee, 2008), while the other
dataset pertains to a study on party-level advantage in U.S. Senate elections
spanning from 1914 to 2010 (Cattaneo et al., 2015). In these datasets, the
forcing variable $x_{i}$ represents the margin of victory of the Democratic
party in a given election, with the outcome variable being the Democratic vote
share in the subsequent election.

Similar to the analysis on simulated datasets, we compare our approach with
the local estimator. The results are reported in Table 5, and the estimated
curves are plotted in Figure 2. For both datasets, the PL method performs
comparably to the local estimator. This finding is consistent with our
observations based on the synthetic datasets, suggesting that PL could serve
as a viable alternative option for RDDs in practice.

Table 5: The performance of the PL method for real data sets.

Method | U.S. House elections | U.S. Senate elections
---|---|---
estimate | s.e. | CL |  $p$-value | estimate | s.e. | CL |  $p$-value
IK1 | 5.808 | 1.083 | 2.124 | 9.9e-8 | 0.082 | 0.008 | 0.016 | 0
PL1 | 7.156 | 1.721 | 3.375 | 3.4e-5 | 0.065 | 0.016 | 0.031 | 4.3e-5
PL2 | 5.493 | 1.167 | 2.289 | 2.8e-6 | 0.074 | 0.013 | 0.026 | 1.4e-8

Figure 2: Estimated treatment effects in the datasets of U.S. House elections
(upper) and U.S. Senate elections (lower).

We lack knowledge of the true treatment effects in the real datasets. To
further assess the performance of our method relative to the local estimator,
we conduct an experiment utilizing the real-world datasets. We set a cutoff at
$c=0.05$ or -0.05 and maintained the original control and treated units,
generating corresponding virtual datasets. Applying our global estimator and
the local estimator to these virtual datasets, we infer the treatment effect
at $c=0.05$ or -0.05 and present the results in Table 6. Our method fails to
reject the null hypothesis that the treatment effect is 0, whereas the local
estimator significantly rejects it. Under the reasonable assumption that the
treatment effect at $c=0.05$ or -0.05 is truly 0, our method offers more
reliable inferences within this context. This investigation highlights the
potential advantage of our method over the local estimator in inference for
RDDs.

Table 6: The performance on the simulated datasets originated from the real
datasets.

cutoff | Method | U.S. House elections | U.S. Senate elections
---|---|---|---
estimate | s.e. | CL |  $p$-value | estimate | s.e. | CL |  $p$-value
0.05 | IK1 | 0.0496 | 0.0076 | 0.0149 | 0.0000 | 2.8820 | 1.1550 | 2.2663 | 0.0127
PL1 | -0.0040 | 0.0142 | 0.0279 | 0.7774 | 1.8499 | 2.0847 | 4.0897 | 0.3750
PL2 | -0.0033 | 0.0146 | 0.0287 | 0.8228 | 1.8639 | 1.3064 | 2.5628 | 0.1539
-0.05 | IK1 | 0.0509 | 0.0072 | 0.0142 | 0.0000 | 2.2451 | 1.1703 | 2.2966 | 0.0553
PL1 | -0.0050 | 0.0149 | 0.0293 | 0.7401 | -0.5513 | 1.9617 | 3.8484 | 0.7787
PL2 | -0.0013 | 0.0147 | 0.0288 | 0.9278 | 1.8973 | 1.2292 | 2.4114 | 0.1229

##  6 Discussion

While local estimators are commonly preferred over global estimators for RDDs,
this paper demonstrates that a semiparametric inference framework provides
efficient global estimation for RDDs. We have applied $P$-spline to represent
the nonparametric part, but other nonparametric methods may be efficient as
well. Further studies are needed. We have examined this inference procedure
within the classical setting of sharp RDDs. It is anticipated that the
semiparametric inference framework can be adapted for use with fuzzy RDDs with
appropriate modifications. Further in-depth analysis is required to explore
this potential extension fully. Another potential extension is the regression
kink design (Card et al., 2015), where Calonico et al. (2014) and Gelman and
Imbens (2019) recommend using local quadratic approach, extending the argument
presented by Hahn et al. (2001). Similar to our study on RDDs, our
semiparametric inference may offer an alternative to local estimators in the
regression kink design. Furthermore, exploring RDDs with covariates (Calonico
et al., 2019) and RDDs with a discrete running variable (Kolesár and Rothe,
2018) could be worthwhile extensions to consider.

## References

  * Bates et al. (2015) Bates, D., Mächler, M., Bolker, B., and Walker, S. (2015), “Fitting Linear Mixed-Effects Models Using lme4,” Journal of Statistical Software, 67, 1–48.
  * Calonico et al. (2019) Calonico, S., Cattaneo, M., Farrell, M., and Titiunik, R. (2019), “Regression Discontinuity Designs using Covariates,” The Review of Economics and Statistics, 101, 442–451.
  * Calonico et al. (2014) Calonico, S., Cattaneo, M. D., and Titiunik, R. (2014), “Robust Nonparametric Confidence Intervals for Regression Discontinuity Designs,” Econometrica, 82, 2295–2326.
  * Calonico et al. (2015) — (2015), “Optimal Data-Driven Regression Discontinuity Plots,” Journal of the American Statistical Association, 110, 1753–1769.
  * Card et al. (2015) Card, D., Lee, D. S., Pei, Z., and Weber, A. (2015), “Inference on Causal Effects in a Generalized Regression Kink Design,” Econometrica, 83, 2453–2483.
  * Cattaneo et al. (2015) Cattaneo, M., Frandsen, B., and Titiunik, R. (2015), “Randomization Inference in the Regression Discontinuity Design: An Application to Party Advantages in the U.S. Senate,” Journal of Causal Inference, 3, 1–24.
  * Cattaneo et al. (2019) Cattaneo, M. D., Idrobo, N., and Titiunik, R. (2019), A Practical Introduction to Regression Discontinuity Designs: Foundations, Cambridge: Cambridge University Press.
  * Claeskens et al. (2009) Claeskens, G., Krivobokova, T., and Opsomer, J. (2009), “Asymptotic properties of penalized spline estimators,” Biometrika, 96, 529–544.
  * Fan (1992) Fan, J. (1992), “Design Adaptive Nonparametric Regression,” Journal of the American Statistical Association, 87, 998–1004.
  * Fan and Gijbels (1996) Fan, J. and Gijbels, I. (1996), Local Polynomial Modelling and Its Applications, London: Chapman & Hall.
  * Gelman and Imbens (2019) Gelman, A. and Imbens, G. W. (2019), “Why High-order Polynomials Should Not Be Used in Regression Discontinuity Designs,” Journal of Business & Economic Statistics, 37, 447–456.
  * Hahn et al. (2001) Hahn, J., Todd, P., and Van der Klaauw, W. (2001), “Identification and Estimation of Treatment Effects With a Regression-Discontinuity Design,” Econometrica, 69, 201–209.
  * Imbens and Kalyanaraman (2012) Imbens, G. and Kalyanaraman, K. (2012), “Optimal Bandwidth Choice for the Regression Discontinuity Estimator,” Review of Economic Studies, 79, 933–959.
  * Imbens and Lemieux (2008) Imbens, G. and Lemieux, T. (2008), “Regression Discontinuity Designs: A Guide to Practice,” Journal of Econometrics, 142, 615–635.
  * Jiang (1998) Jiang, J. (1998), “Asymptotic properties of the empirical BLUP and BLUE in mixed linear models,” Statistica Sinica, 8, 861–885.
  * Kolesár and Rothe (2018) Kolesár, M. and Rothe, C. (2018), “Inference in Regression Discontinuity Designs with a Discrete Running Variable,” American Economic Review, 108, 2277–2304.
  * Lee (2008) Lee, D. (2008), “Randomized Experiments from Non-Random Selection in U.S. House Elections,” Journal of Econometrics, 142, 675–697.
  * Lee and Lemieux (2010) Lee, D. and Lemieux, T. (2010), “Regression Discontinuity Designs in Economics,” Journal of Economic Literature, 48, 281–355.
  * Miller (1977) Miller, J. (1977), “Asymptotic properties of maximum likelihood estimates in the mixed model of the analysis of variance,” The Annals of Statistics, 5, 746–762.
  * Patterson and Thompson (1971) Patterson, H. D. and Thompson, R. (1971), “Recovery of inter-block information when block sizes are unequal,” Biometrika, 58, 545–554.
  * Pinheiro et al. (2023) Pinheiro, J., Bates, D., and R Core Team (2023), nlme: Linear and Nonlinear Mixed Effects Models, r package version 3.1-164.
  * Ruppert (2002) Ruppert, D. (2002), “Selecting the number of knots for penalized splines,” Journal of Computational and Graphical Statistics, 11, 735–757.
  * Ruppert et al. (2003) Ruppert, D., Wand, M. P., and Carroll, R. J. (2003), Semiparametric Regression, Cambridge: Cambridge University Press.
  * Thistlewaite and Campbell (1960) Thistlewaite, D. and Campbell, D. (1960), “Regression-Discontinuity Analysis: An Alternative to the Ex-Post Facto Experiment,” Journal of Educational Psychology, 51, 309–317.

## Appendix A: The proof of Theorem 2

Let $\nu_{1}=n-K_{r}$, $\nu_{2}=K_{r}$, where $K_{r}=\text{rank}(\mathbf{Z})$.
We define the matrices $\mathbf{G}_{1}=\sigma^{2}\mathbf{I}_{n}$ and
$\mathbf{G}_{2}=\sigma_{\gamma}^{2}\mathbf{Z}\mathbf{Z}^{\top}$, and introduce
the $2\times 2$ matrix $\mathbf{C}$ with its $(i,j)$-element defined as
follows:

$[\mathbf{C}]_{ij}=\lim_{n\rightarrow\infty}\frac{\text{tr}[\mathbf{V}^{-1}\mathbf{G}_{i}\mathbf{V}^{-1}\mathbf{G}_{j}]}{\sqrt{\nu_{i}\nu_{j}}}.$

Denote
$\boldsymbol{\eta}=(\sigma^{2},\sigma_{\gamma}^{2},\boldsymbol{\theta}^{\top})^{\top}$
and
$\hat{\boldsymbol{\eta}}=(\hat{\sigma}^{2},\hat{\sigma}_{\gamma}^{2},\hat{\boldsymbol{\theta}}^{\top})^{\top}$.
Miller (1977) established the asymptotic properties of maximum likelihood
estimates $\hat{\boldsymbol{\eta}}$, in a general mixed model under the
following conditions:
A1: $K\rightarrow\infty$ as $n\rightarrow\infty$,
A2: $\lim_{n\rightarrow\infty}K/n=0$,
A3: $\lim_{n\rightarrow\infty}K_{r}/K>0$,
A4: $\mathbf{J}$ is positive definite,
A5: $\mathbf{C}$ is positive definite.
Theorem 3.1 of Miller (1977) demonstrates that under conditions A1-A5, as
$n\rightarrow\infty$,

$\sqrt{n}(\hat{\boldsymbol{\eta}}-\boldsymbol{\eta})\rightarrow_{d}\mathcal{N}(0,\text{diag}\\{\mathbf{C}^{-1},\mathbf{J}^{-1}\\}).$

Theorem 2 is an application of Theorem 3.1 of Miller (1977). It is sufficient
to verify that the conditions (A1), (A2), and (A3) as stated in Theorem 2 are
satisfied. It remains to demonstrate the fulfillment of (A4) and (A5).

First, we analyze the $(p+2)\times(p+2)$ matrix

$[\mathbf{J}]_{ij}=\lim_{n\rightarrow\infty}\mathbf{U}^{\top}\mathbf{V}^{-1}\mathbf{U}/n.$

Simple derivation shows

$\mathbf{V}^{-1}=\sigma^{-2}\mathbf{I}_{n}-\sigma^{-2}\mathbf{Z}(\sigma_{\gamma}^{-2}\sigma^{2}\mathbf{I}_{K}+\mathbf{Z}^{\top}\mathbf{Z})^{-1}\mathbf{Z}^{\top}.$

Note
$\mathbf{U}^{\top}\mathbf{Z}=\left(\begin{array}[]{c}\mathbf{w}^{\top}\mathbf{Z}\\\
\mathbf{X}^{\top}\mathbf{Z}\end{array}\right)$, and
$\mathbf{U}^{\top}\mathbf{U}=\left(\begin{array}[]{cc}\mathbf{w}^{\top}\mathbf{w}&\mathbf{w}^{\top}\mathbf{X}\\\
\mathbf{X}^{\top}\mathbf{w}&\mathbf{X}^{\top}\mathbf{X}\end{array}\right)$. We
have that

$\displaystyle\sigma^{2}\mathbf{U}^{\top}\mathbf{V}^{-1}\mathbf{U}=$
$\displaystyle\mathbf{U}^{\top}\mathbf{U}-\mathbf{U}^{\top}\mathbf{Z}(\sigma_{\gamma}^{-2}\sigma^{2}\mathbf{I}_{K}+\mathbf{Z}^{\top}\mathbf{Z})^{-1}\mathbf{Z}^{\top}\mathbf{U}$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\mathbf{w}^{\top}\mathbf{w}&\mathbf{w}^{\top}\mathbf{X}\\\
\mathbf{X}^{\top}\mathbf{w}&\mathbf{X}^{\top}\mathbf{X}\end{array}\right)-\left(\begin{array}[]{c}\mathbf{w}^{\top}\mathbf{Z}\\\
\mathbf{X}^{\top}\mathbf{Z}\end{array}\right)\left(\frac{\sigma^{2}}{\sigma_{\gamma}^{2}}\mathbf{I}_{K}+\mathbf{Z}^{\top}\mathbf{Z}\right)^{-1}(\mathbf{Z}^{\top}\mathbf{w},\mathbf{Z}^{\top}\mathbf{X}).$
(23)

Denoting
$\boldsymbol{\mu}_{x}=\lim_{n\rightarrow\infty}n^{-1}\mathbf{X}^{\top}\mathbf{1}$,
$\boldsymbol{\Sigma}_{x}=\lim_{n\rightarrow\infty}n^{-1}\mathbf{X}^{\top}\mathbf{X}$,
$\boldsymbol{\Sigma}_{xz}=\lim_{n\rightarrow\infty}n^{-1}\mathbf{X}^{\top}\mathbf{Z}$.
and noting
$\lim_{n\rightarrow\infty}n^{-1}\mathbf{w}^{\top}\mathbf{w}=\lim_{n\rightarrow\infty}n_{1}/n=\rho_{1}$,
$\lim_{n\rightarrow\infty}n^{-1}\mathbf{w}^{\top}\mathbf{X}=\rho_{1}\boldsymbol{\mu}_{x}^{\top}$,
and
$\lim_{n\rightarrow\infty}n^{-1}\mathbf{w}^{\top}\mathbf{Z}=\rho_{1}\boldsymbol{\mu}_{z}^{\top}$,
(Appendix A: The proof of Theorem 2) follows that

$\displaystyle\lim_{n\rightarrow\infty}n^{-1}\mathbf{U}^{\top}\mathbf{V}^{-1}\mathbf{U}=$
$\displaystyle\sigma^{-2}\left[\left(\begin{array}[]{cc}\rho_{1}&\rho_{1}\boldsymbol{\mu}_{x}^{\top}\\\
\rho_{1}\boldsymbol{\mu}_{x}&\boldsymbol{\Sigma}_{x}\end{array}\right)-\left(\begin{array}[]{c}\rho_{1}\boldsymbol{\mu}_{z}^{\top}\\\
\boldsymbol{\Sigma}_{xz}\end{array}\right)\boldsymbol{\Sigma}_{z}^{-1}\left(\rho_{1}\boldsymbol{\mu}_{z},\boldsymbol{\Sigma}_{xz}^{\top}\right)\right]$
(28) $\displaystyle=$
$\displaystyle\sigma^{-2}\left(\begin{array}[]{cc}\rho_{1}-\rho_{1}^{2}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}&\rho_{1}\boldsymbol{\mu}_{x}^{\top}-\rho_{1}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}\\\
\rho_{1}\boldsymbol{\mu}_{x}-\rho_{1}\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}&\boldsymbol{\Sigma}_{x}-\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}\end{array}\right).$
(31)

For any vector $\mathbf{u}=(u_{1},\mathbf{u}_{2}^{\top})^{\top}\neq
0\in\mathbb{R}^{p+2}$, we have that

$\displaystyle\mathbf{u}^{\top}\left(\begin{array}[]{cc}\rho_{1}-\rho_{1}^{2}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}&\rho_{1}\boldsymbol{\mu}_{x}^{\top}-\rho_{1}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}\\\
\rho_{1}\boldsymbol{\mu}_{x}-\rho_{1}\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}&\boldsymbol{\Sigma}_{x}-\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}\end{array}\right)\mathbf{u}$
(34) $\displaystyle=$ $\displaystyle
u_{1}^{2}(\rho_{1}-\rho_{1}^{2}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z})+\mathbf{u}_{2}^{\top}\left(\boldsymbol{\Sigma}_{x}-\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}\right)\mathbf{u}_{2}+2u_{1}\rho_{1}\left(\boldsymbol{\mu}_{x}^{\top}-\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}\right)\mathbf{u}_{2}$
$\displaystyle=$
$\displaystyle\left[\rho_{1}u_{1}+\left(\boldsymbol{\mu}_{x}^{\top}-\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}\right)\mathbf{u}_{2}\right]^{2}+u_{1}^{2}(\rho_{1}-\rho_{1}^{2}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z})-\rho_{1}^{2}u_{1}^{2}$
$\displaystyle+\mathbf{u}_{2}^{\top}\left[\boldsymbol{\Sigma}_{x}-\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}-(\boldsymbol{\mu}_{x}-\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z})(\boldsymbol{\mu}_{x}^{\top}-\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top})\right]\mathbf{u}_{2}.$
(35)

In the right-hand of (34),
$\left[\rho_{1}u_{1}+\left(\boldsymbol{\mu}_{x}^{\top}-\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}\right)\mathbf{u}_{2}\right]^{2}\geq
0$ and

$\mathbf{u}_{2}^{\top}\left[\boldsymbol{\Sigma}_{x}-\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}-(\boldsymbol{\mu}_{x}-\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z})(\boldsymbol{\mu}_{x}^{\top}-\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top})\right]\mathbf{u}_{2}\geq
0$

due to
$\boldsymbol{\Sigma}_{x}-\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}>(\boldsymbol{\mu}_{x}-\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z})(\boldsymbol{\mu}_{x}^{\top}-\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top})$,
a result demonstrated in Lemma 1 below. Here, for matrices $\mathbf{A}$ and
$\mathbf{B}$, $\mathbf{A}>\mathbf{B}$ means that the matrix
$\mathbf{A}-\mathbf{B}$ is positive. To ensure that
$\lim_{n\rightarrow\infty}n^{-1}\mathbf{U}^{\top}\mathbf{V}^{-1}\mathbf{U}$ is
positive, we only need the following inequality to hold:

$\displaystyle\rho_{1}-\rho_{1}^{2}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}>$
$\displaystyle\rho_{1}^{2}$ (36)

Clearly, (36) is satisfied under the part (d) of Assumption 4.

Next, Let us check Condition (A5). Denote
$\mathbf{P}_{z}=\mathbf{Z}(\sigma_{\gamma}^{-2}\sigma^{2}\mathbf{I}_{K}+\mathbf{Z}^{\top}\mathbf{Z})^{-1}\mathbf{Z}^{\top}$,
implying $\mathbf{V}^{-1}=\sigma^{-2}(\mathbf{I}_{n}-\mathbf{P}_{z})$. Thus,
we have that

$\displaystyle\mathbf{V}^{-2}=$
$\displaystyle\sigma^{-4}\left(\mathbf{I}_{n}-2\mathbf{P}_{z}+\mathbf{P}_{z}^{2}\right),$
$\displaystyle\mathbf{Z}^{\top}\mathbf{V}^{-1}\mathbf{Z}=$
$\displaystyle\sigma^{-2}\left(\mathbf{Z}^{\top}\mathbf{Z}-\mathbf{Z}^{\top}\mathbf{P}_{z}\mathbf{Z}\right)=\sigma^{-2}\sigma_{\gamma}^{-2}(\sigma_{\gamma}^{-2}\sigma^{2}\mathbf{I}_{K}+\mathbf{Z}^{\top}\mathbf{Z})^{-1}\mathbf{Z}^{\top}\mathbf{Z},$
$\displaystyle\mathbf{Z}\mathbf{V}^{-2}\mathbf{Z}^{\top}=$
$\displaystyle\sigma^{-4}\left(\mathbf{Z}^{\top}\mathbf{Z}-2\mathbf{Z}^{\top}\mathbf{P}_{z}\mathbf{Z}+\mathbf{Z}^{\top}\mathbf{P}_{z}^{2}\mathbf{Z}\right).$

Based of the equations above, the elements of $\mathbf{C}$ are represented by

$\displaystyle[\mathbf{C}]_{11}=$
$\displaystyle\lim_{n\rightarrow\infty}\frac{\sigma^{4}}{n-K_{r}}\text{tr}(\mathbf{V}^{-2})$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{n-K_{r}}\left(n-2\text{tr}\left[\mathbf{P}_{z}\right]+\text{tr}\left[\mathbf{P}_{z}^{2}\right]\right)=\lim_{n\rightarrow\infty}\frac{n-K}{n-K_{r}}=1,$
(37) $\displaystyle[\mathbf{C}]_{22}=$
$\displaystyle\lim_{n\rightarrow\infty}\frac{\sigma_{\gamma}^{4}}{K_{r}}\text{tr}\left[(\mathbf{Z}\mathbf{V}^{-1}\mathbf{Z}^{\top})^{2}\right]$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{K_{r}}\text{tr}\left[\left((\sigma_{\gamma}^{-2}\sigma^{2}\mathbf{I}_{K}+\mathbf{Z}^{\top}\mathbf{Z})^{-1}\mathbf{Z}^{\top}\mathbf{Z}\right)^{2}\right]=\lim_{n\rightarrow\infty}\frac{K}{K_{r}}=1/r,$
(38) $\displaystyle[\mathbf{C}]_{12}=$
$\displaystyle\lim_{n\rightarrow\infty}\frac{\sigma^{2}\sigma_{\gamma}^{2}}{\sqrt{K_{r}(n-K_{r})}}\text{tr}(\mathbf{Z}\mathbf{V}^{-2}\mathbf{Z}^{\top})$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\frac{\sigma^{2}\sigma_{\gamma}^{-2}\text{tr}\left[\mathbf{Z}^{\top}\mathbf{Z}-2\mathbf{Z}^{\top}\mathbf{P}_{z}\mathbf{Z}+\mathbf{Z}^{\top}\mathbf{P}_{z}^{2}\mathbf{Z}\right]}{\sqrt{K_{r}(n-K_{r})}}=\lim_{n\rightarrow\infty}\frac{K}{\sqrt{K_{r}(n-K_{r})}}=0,$
(39)

where the last step of (37) is from Part (b) of Assumption 4, the last step of
(38) is from Part (c) of Assumption 4, and the last step of (39) is from Parts
(b) & (c) of Assumption 4. Therefore, it establishes A5 as desired.

######  Lemma 1.

$\boldsymbol{\Sigma}_{x}-\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}>(\boldsymbol{\mu}_{x}-\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z})(\boldsymbol{\mu}_{x}^{\top}-\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top})$

###### Proof.

The inequality is equivalent to

$\boldsymbol{\Sigma}_{x}-\boldsymbol{\mu}_{x}\boldsymbol{\mu}_{x}^{\top}+\boldsymbol{\mu}_{x}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}+\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}\boldsymbol{\mu}_{x}^{\top}-\boldsymbol{\Sigma}_{xz}(\boldsymbol{\Sigma}_{z}^{-1}+\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1})\boldsymbol{\Sigma}_{xz}^{\top}>0.$

Because
$\boldsymbol{\Sigma}_{x}-\boldsymbol{\mu}_{x}\boldsymbol{\mu}_{x}^{\top}\geq
0$, it is necessary to prove that
$\boldsymbol{\mu}_{x}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}+\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}\boldsymbol{\mu}_{x}^{\top}-\boldsymbol{\Sigma}_{xz}(\boldsymbol{\Sigma}_{z}^{-1}+\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1})\boldsymbol{\Sigma}_{xz}^{\top}>0$.
This is equivalent to checking

$\boldsymbol{v}^{\top}[\boldsymbol{\mu}_{x}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\Sigma}_{xz}^{\top}+\boldsymbol{\Sigma}_{xz}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}\boldsymbol{\mu}_{x}^{\top}-\boldsymbol{\Sigma}_{xz}(\boldsymbol{\Sigma}_{z}^{-1}+\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1})\boldsymbol{\Sigma}_{xz}^{\top}]\boldsymbol{v}\geq
0.$ (40)

Define

$f(\mathbf{x})=\boldsymbol{v}^{\top}\boldsymbol{\mu}_{x}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\mathbf{x}^{\top}+\mathbf{x}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}\boldsymbol{\mu}_{x}^{\top}\boldsymbol{v}-\mathbf{x}(\boldsymbol{\Sigma}_{z}^{-1}+\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1})\mathbf{x}^{\top}.$

The inequality (40) becomes
${f}(\boldsymbol{v}^{\top}\boldsymbol{\Sigma}_{xz})\geq 0$. Setting the first
derivative of ${f}(\mathbf{x})$ to 0, we need to solve the following equation:

$(\boldsymbol{\Sigma}_{z}^{-1}+\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1})\mathbf{x}^{\top}=\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}\boldsymbol{\mu}_{x}^{\top}\boldsymbol{v}.$
(41)

Let
$\mathbf{A}=\boldsymbol{\Sigma}_{z}^{-1}+\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}$.
Since
$\mathbf{I}+\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\boldsymbol{\mu}_{z}>0$,
$\mathbf{A}$ is invertible by the Sherman-Morrison-Woodbury formula. Thus, we
obtain the unique solution to (41) as follows:

$\tilde{\mathbf{x}}=\boldsymbol{v}^{\top}\boldsymbol{\mu}_{x}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma}_{z}^{-1}\mathbf{A}^{-1}.$

Moreover,

${f}(\tilde{\mathbf{x}})=(\boldsymbol{v}^{\top}\boldsymbol{\mu_{x}}\boldsymbol{\mu}_{z}^{\top}\boldsymbol{\Sigma_{z}}^{-1})\mathbf{A}^{-1}(\boldsymbol{\Sigma_{z}}^{-1}\boldsymbol{\mu_{z}}\boldsymbol{\mu}_{x}^{\top}\boldsymbol{v})>0.$

Theorefore, the lemma is proved. ∎