text:
\section{Introduction}

\subsection{Preamble}
\label{pre}

Panel data analysis is an important topic in statistics and econometrics.  The
traditional approach for analysing panel data assumes all individuals share the
same unknown parameters, and uses cluster effects to account for the difference
between individuals.  For example, when the linear models are
used, the coefficients of the covariates are assumed to be the same across all
individuals, i.e.
$$
y_{it} = {\bf X}_{it}^{{\mbox{\rm\tiny T}}} \mbox{\boldmath $\beta$} + \epsilon_{it},
\quad
i=1, \ \cdots, \ m;  \  \  t=1, \ \cdots, \ T,
$$
where $y_{it}$ and ${\bf X}_{it}$, a $(p+1)$-dimensional vector, are respectively
the $t$th observations of the response variable and covariate of the $i$th
individual.  $\epsilon_{it}$, $t=1, \ \ldots, \ T$, are correlated for any
given $i$, and the cluster effects are included in $\epsilon_{it}$.  See
\cite{hsiao14} and the reference therein.  Whilst this modelling idea is
useful when the global trend of the impact of a covariate on the response
variable is of our main interest, it does not tell us anything about the
individual attributes which are sometimes very important.

In order to explore the individual attributes, we need to make them more
concrete and distinctive in modelling.  A simple approach to do so would be
using
\begin{equation}
y_{it} = {\bf X}_{it}^{{\mbox{\rm\tiny T}}} \mbox{\boldmath $\beta$}_i + \epsilon_{it},
\quad
i=1, \ \cdots, \ m;  \  \  t=1, \ \cdots, \ T,
\label{eq1}
\end{equation}
to fit the data.  However, this modelling approach would result in $m(p+1)$
unknown coefficients to estimate, which is too many, because $m$ is
usually of the magnitude of hundreds, or even more, in practice.  This
modelling also ignores the similarity which may exist among some individuals.
Such similarity may have very important practical meaning, and could lead to
some important findings in practice.  In addition to that, statistically
speaking, the modelling, like (\ref{eq1}) without any conditions imposed, would
also pay a price on variance side of the estimators resulted because the
available information is not used up.

In order to explore the individual attributes and account for the similarity
among some individuals at the same time, \cite{ke15} proposed a penalised
likelihood/least squares based approach to pursue the homogeneity in the
linear models, i.e. (\ref{eq1}), used for panel data analysis, under the
framework of treating homogeneity as a kind of sparsity.   Regression under
homogeneity condition has also been studied by quite a few recent works, e.g.
 \cite{tibshirani05,friedman07,bondell2008simultaneous,jiangqian13}
, and the references therein.  Like \cite{ke15},
the methods in these works are all based on penalised likelihood/least
squares.  \cite{ke16} took a different approach, they formulated the
homogeneity pursuit problem as a problem of change point detection and applied
the binary segmentation approach to identify the homogeneity in the
linear models with interactive effects.

The existing literature about homogeneity pursuit mainly focuses on the linear
models.  It is well known that the linearity condition may not hold for many
datasets, and the exploration of linear relationship is not sufficient in many
cases.  As a consequence, the semiparametric modelling is becoming more and
more useful in panel data analysis.  Among various semiparametric models, the
single index models
have many advantages, and are a very successful tool in data analysis, see
 \cite{hardle1989investigating,carroll97,yuruppert02,zhuxue06,xia08,penghuang11,zhumiaopeng12,guo2016dynamic},
 and the
reference therein.  In this paper, we are going to investigate the homogeneity
pursuit in the single index models used for panel data analysis.  The detailed
definition of the models we are going to address in this paper is given in
Section \ref{mo0}

\subsection{The single index models with homogeneity structure}
\label{mo0}

Let $y_{it}$ and ${\bf X}_{it}$, a $(p+1)$-dimensional vector, be
respectively the $t$th observations of the response variable and covariate of
the $i$th individual, $i=1, \ \cdots, \ m; \ \  t=1, \ \cdots, \ T$. We
consider the models
\begin{equation}
y_{it}
=
g_i({\bf X}_{it}^{\rm T}\mbox{\boldmath $\beta$}_i) + \epsilon_{it},
\quad
i=1, \ \cdots, \ m; \ \  t=1, \ \cdots, \ T,
\label{model}
\end{equation}
where
\begin{equation}
g_i(\cdot)
=
\left\{
\begin{array}{ll}
\ g_{(1)}(\cdot) & \mbox{ when }  i \in G_{1,1},
\\
\ g_{(2)}(\cdot) & \mbox{ when }  i \in G_{1,2},
\\
\quad \vdots & \quad \vdots
\\
\ g_{(H_1)}(\cdot) & \mbox{ when }  i \in G_{1,H_1},
\end{array}
\right.
\quad
\beta_{ij}
=
\left\{
\begin{array}{ll}
\ \beta_{(1)} & \mbox{ when }  (i, j) \in G_{2,1},
\\
\ \beta_{(2)} & \mbox{ when }  (i, j) \in G_{2,2},
\\
\quad \vdots & \quad \vdots
\\
\ \beta_{(H_2)} & \mbox{ when }  (i, j) \in G_{2,H_2},
\end{array}
\right.
\label{con}
\end{equation}
$\mathbb{G}_1=\{G_{1,k}: \ k = 1, \ \cdots, \ H_1\}$ is a partition of set
$\{1, \ \cdots, \ m\}$, $\mathbb{G}_2=\{G_{2,k}: \ k = 1, \ \cdots, \ H_2\}$ is a partition
of set $\{(i,j): \ i = 1, \ \cdots, \ m; \ j=1, \ \cdots, \ p\}$, $\beta_{ij}$
is the $(j+1)$th component of $\mbox{\boldmath $\beta$}_i$,
and
$$
E(\epsilon_{it}| {\bf X}_{it}) = 0,
\quad
\mbox{var}(\epsilon_{it}| {\bf X}_{it}) = \sigma^2.
$$
The condition (\ref{con}) is the homogeneity structure of the standard single
index models for panel data analysis.  $\{G_{1,k}: \ k = 1, \ \cdots, \ H_1\}$
and $\{G_{2,k}: \ k = 1, \ \cdots, \ H_2\}$ are unknown partitions.  $H_1$ and
$H_2$ are unknown integers, $H_1$ is much smaller than $m$, $H_2$ is much
smaller than $mp$.  $g_{(k)}(\cdot)$, $k=1, \ \cdots, \ H_1$, are unknown
functions to be estimated, and $\beta_{(k)}$, $k=1, \ \cdots, \ H_2$, are
unknown parameters to be estimated.

Let $\beta_{i0}=1$ be the first component of $\mbox{\boldmath $\beta$}_i$.  In the
literature, the most commonly used identification condition for the single
index models is $\|\mbox{\boldmath $\beta$}_i\|=1$ and $\beta_{i0}>0$, or $\beta_{i0}=1$.  We
choose the latter in this paper.

The models (\ref{model}) together with (\ref{con}) show that the homogeneity
pursuit in the single index models for panel data analysis is even more
important than that in the linear models, this is because we would have to
estimate $m$ unknown functions and $mp$ unknown parameters in order to explore
the individual attributes, if the homogeneity pursuit is not conducted.
However, if the homogeneity pursuit is conducted, we only need to estimate
$H_1$, much smaller than $m$, unknown functions and $H_2$, much smaller than
$mp$, unknown parameters when the homogeneity exists.  Even without taking into
account the benefit resulted from the homogeneity pursuit for the parametric
part of the models, just for the part of unknown functions alone, to estimate
much fewer functions would make a big difference in the obtained estimators, in
terms of the stability of the estimators.

The rest of the paper is organized as follows.  We begin in Section \ref{est}
with a description of the proposed estimation procedure which is embedded with
a binary segmentation based homogeneity pursuit.  The asymptotic properties of
the proposed estimators are presented in Section \ref{asy}.  The performance of
the proposed estimation procedure and homogeneity pursuit method, when sample
size is finite, are assessed by simulation studies in Section \ref{simu}.  In
Section \ref{re}, applying the single index models (\ref{model}) together with
the homogeneity structure (\ref{con}) to the 49 Industry Portfolios data set,
which can be freely downloaded from Kenneth French's website
\begin{center}
{\tt
http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data\_library.html},
\end{center}
and the UK climate data, which can be freely downloaded from
\begin{center}
{\tt http://www.metoffice.gov.uk/public/weather/climate-historic},
\end{center}
we will show the advantages of the proposed statistical methodology.  We leave
all technical proofs of the asymptotic properties in the Appendix.

\section{Estimation procedure}
\label{est}

\subsection{Estimation method}
\label{est1}

Our approach to deal with the unknown functions $g_i(\cdot)$,
$i=1, \ \cdots, \ m$, in (\ref{model}) is based on the B-Spline.  To achieve
the best result
for the homogeneity pursuit, we have to decompose all $g_i(\cdot)$s by the same
B-Spline basis, ${\bf B}(\cdot) = (B_1(\cdot), \ \cdots, \ B_K(\cdot))^{{\mbox{\rm\tiny T}}}$.

For each $i$, $i = 1, \ \cdots, \ m$, let $\tilde{\mbox{\boldmath $\beta$}}_i$ be the estimate
of $\mbox{\boldmath $\beta$}_i$ obtained, based on the observations for the $i$th individual,
by a standard estimation procedure for the single index models, e.g. the
method in \cite{yuruppert02} or in \cite{hardle1989investigating}, and
$$
a
=
\min\limits_{1 \leq i \leq m} \min\limits_{1 \leq t \leq T}
{\bf X}_{it}^{{\mbox{\rm\tiny T}}} \tilde{\mbox{\boldmath $\beta$}}_i,
\quad
b
=
\max\limits_{1 \leq i \leq m} \max\limits_{1 \leq t \leq T}
{\bf X}_{it}^{{\mbox{\rm\tiny T}}} \tilde{\mbox{\boldmath $\beta$}}_i.
$$
We use the B-Spline basis of order $s$ in this paper, and the basis,
${\bf B}(\cdot)$, is formed by the equally spaced knots, $\tau_k$,
$k=0, \ \cdots, \ K-s+1$, on the interval $[a, \ b]$, with $\tau_0 = a$ and
$\tau_{K-s+1} = b$.  Based on the basis ${\bf B}(\cdot)$, $g_i(\cdot)$ can be
decomposed as
\begin{equation}
g_i(\cdot) \approx {\bf B}(\cdot)^{{\mbox{\rm\tiny T}}} \mbox{\boldmath $\theta$}_i,
\label{app}
\end{equation}
where $\mbox{\boldmath $\theta$}_i = (\theta_{i1}, \ \cdots, \ \theta_{iK})^{{\mbox{\rm\tiny T}}}$.  So, to get
the estimator of $g_i(\cdot)$, we only need to get the estimator of $\mbox{\boldmath $\theta$}_i$.

Our estimation procedure for $\mbox{\boldmath $\theta$}_i$ and $\mbox{\boldmath $\beta$}_i$, $i=1, \ \cdots, \ m$,
consists of three stages: in the first stage, for each $i$, we estimate
$\mbox{\boldmath $\theta$}_i$ and $\mbox{\boldmath $\beta$}_i$ only based
on the observations for the $i$th individual, and treat the obtained
estimators as initial estimators; we identify, in the second stage, the
homogeneity structure in the $\mbox{\boldmath $\theta$}_i$s and $\mbox{\boldmath $\beta$}_i$s based on the initial
estimators obtained in the first stage; in the final stage, we estimate the
$\mbox{\boldmath $\theta$}_i$s and $\mbox{\boldmath $\beta$}_i$s under the identified homogeneity structure.

We now present the details of the estimation procedure.
\begin{enumerate}
\item[{\bf Stage 1}] ({\it Initial Estimation}).  Let
$\overline\mbox{\boldmath $\beta$}_i=(\beta_{i1}, \ \cdots, \ \beta_{ip})^{\rm T}$, which
is $\mbox{\boldmath $\beta$}_i$ with the first component, which is always $1$, being dropped.
For each $i$, based on the observations for the $i$th individual, approximating
$g_i(\cdot)$ by its decomposition (\ref{app}) and applying the least squares
estimation method, we have the following objective function
\begin{equation}
\sum_{t=1}^T
\Big(
y_{it} - {\bf B}^{\rm T}({\bf X}_{it}^{\rm T}\mbox{\boldmath $\beta$}_i)\mbox{\boldmath $\theta$}_i
\Big)^2.
\label{eqn:min1}
\end{equation}
Minimise (\ref{eqn:min1}) with respect to
$(\overline\mbox{\boldmath $\beta$}_i^{{\mbox{\rm\tiny T}}}, \ \mbox{\boldmath $\theta$}_i^{{\mbox{\rm\tiny T}}})$, and denote the resulting minimiser by
$(\tilde{\mbox{\boldmath $\beta$}}_i^{{\mbox{\rm\tiny T}}}, \ \tilde{\mbox{\boldmath $\theta$}}_i^{{\mbox{\rm\tiny T}}})$.  We will show how to
conduct the minimisation in Section \ref{est2}.

\item[{\bf Stage 2}] ({\it Homogeneity Pursuit}).  Let $\tilde{\beta}_{ij}$ be
the $j$th component of
$\tilde{\mbox{\boldmath $\beta$}}_i$, we sort $\tilde{\beta}_{ij}$, $i=1, \ \cdots, \ m$,
$j=1, \ \cdots, \ p$, in ascending order, and denote them by
$$
b_{(1)}
\leq
\cdots
\leq
b_{(mp)}.
$$
We use $R_{ij}$ to denote the rank of $\tilde{\beta}_{ij}$.  Identifying the
homogeneity among $\tilde{\beta}_{ij}$, $i=1, \ \cdots, \ m$,
$j=1, \ \cdots, \ p$, is equivalent to detecting the change points among
$b_{(l)}$, $l=1, \ \cdots, \ mp$.  To this end, we apply the Binary
Segmentation algorithm as follows.

For any $1 \leq i <  j \leq mp$, let
$$
\Delta_{ij}(\kappa)
=
\sqrt{\frac{(j-\kappa)(\kappa-i+1)}{j-i+1}}
\left|
\frac{\sum_{l=\kappa+1}^j b_{(l)}}{j-\kappa}
-
\frac{\sum_{l=i}^{\kappa} b_{(l)}}{\kappa-i+1}
\right|.
$$
Given a threshold $\delta$, the Binary Segmentation algorithm to detect the
change points works as follows
\begin{enumerate}
\item[(1)] Find $\hat{k}_1$ such that
$$
\Delta_{1,mp}(\hat{k}_1)
=
\max\limits_{1 \leq \kappa < mp} \Delta_{1,mp}(\kappa).
$$
If $\Delta_{1,mp}(\hat{k}_1) \leq \delta$, there is no change point among
$b_{(l)}$, $l=1, \ \cdots, \ mp$, and the process of detection ends.
Otherwise, add $\hat{k}_1$ to the set of change points and divide the region
$\{\kappa:  \ 1 \leq \kappa \leq mp \}$ into two subregions:
$\{\kappa:  \ 1 \leq \kappa \leq \hat{k}_1 \}$ and
$\{\kappa:  \ \hat{k}_1+1 \leq \kappa \leq mp \}$.

\smallskip
\item[(2)] Detect the change points in the two subregions obtained in (1),
respectively.  Let us deal with the region $\{\kappa:  \ 1 \leq \kappa \leq
\hat{k}_1 \}$ first.  Find $\hat{k}_2$ such that
$$
\Delta_{1,\hat{k}_1}(\hat{k}_2)
=
\max\limits_{1 \leq \kappa < \hat{k}_1} \Delta_{1,\hat{k}_1}(\kappa).
$$
If $\Delta_{1,\hat{k}_1}(\hat{k}_2) \leq \delta$, there is no change point in
the region $\{\kappa:  \ 1 \leq \kappa \leq \hat{k}_1 \}$.  Otherwise,
add $\hat{k}_2$ to the set of change points and divide the region
$\{\kappa:  \ 1 \leq \kappa \leq \hat{k}_1 \}$ into two subregions:
$\{\kappa:  \ 1 \leq \kappa \leq \hat{k}_2 \}$ and
$\{\kappa:  \ \hat{k}_2+1 \leq \kappa \leq \hat{k}_1 \}$.  For the region
$\{\kappa:  \ \hat{k}_1+1 \leq \kappa \leq mp \}$, we find $\hat{k}_3$ such
that
$$
\Delta_{\hat{k}_1+1,mp}(\hat{k}_3)
=
\max\limits_{\hat{k}_1+1 \leq \kappa < mp} \Delta_{\hat{k}_1+1,mp}(\kappa).
$$
If $\Delta_{\hat{k}_1+1,mp}(\hat{k}_3) \leq \delta$, there is no change point in the region
$\{\kappa:  \ \hat{k}_1+1 \leq \kappa \leq mp \}$.  Otherwise, add $\hat{k}_3$
to the set of change points and divide the region $\{\kappa:  \ \hat{k}_1+1
\leq \kappa \leq mp \}$ into two subregions:
$\{\kappa:  \ \hat{k}_1+1 \leq \kappa \leq \hat{k}_3 \}$ and
$\{\kappa:  \ \hat{k}_3+1 \leq \kappa \leq mp \}$.
\smallskip
\item[(3)] For each subregion obtained in (2), we do exactly the same as that
for the subregion $\{\kappa:  \ 1 \leq \kappa \leq \hat{k}_1 \}$ or
$\{\kappa:  \ \hat{k}_1+1 \leq \kappa \leq mp \}$ in (2), and
keep doing so until there is no subregion containing any change point.
\end{enumerate}

We sort the estimated change point locations in ascending order and denote them
by
$$
\hat{k}_{(1)} < \hat{k}_{(2)} < \cdots < \hat{k}_{(\hat{H}_{-1})},
$$
where $\hat{H}_{-1}$ is the number of change points detected.  In addition, we
denote $\hat{k}_{(0)} = 0$, $\hat{H}_2 = \hat{H}_{-1} + 1$, and
$\hat{k}_{(\hat{H}_2)} = mp$.

We use $\hat{H}_2$ to estimate $H_2$.  Let
$$
\hat{G}_{2, s}
=
\{(i, j): \ \hat{k}_{(s-1)} < R_{ij} \leq \hat{k}_{(s)}\},
\quad
1 \leq s \leq \hat{H}_2,
$$
we use $\left\{\hat{G}_{2, s}: \ 1 \leq s \leq \hat{H}_2 \right\}$
to estimate the partition $\{G_{2,s}: \ 1 \leq s \leq H_2\}$.  We consider all
the $\beta_{ij}$s with the subscript $(i,j)$ in the same member of the
estimated partition having the same value.

Let $\tilde{\theta}_{ij}$ be the $j$th component of $\tilde{\mbox{\boldmath $\theta$}}_i$.  Doing
exactly the same to $\tilde{\theta}_{ij}$, $i=1, \ \cdots, \ m$, $j=1, \
\cdots, \ K$, we get a partition $\{\hat G_{1,1}, \ \cdots, \ \hat G_{1,\hat{H}_1}\}$ of
$\{(i, j): \ i=1, \ \cdots, \ m; \ j=1, \ \cdots, \ K\}$.  We consider all the
$\theta_{ij}$s with subscript $(i,j)$ in the same member of the estimated
partition having the same value.

\item[{\bf Stage 3}] ({\it Final Estimation}).  Let
$L(\eta_1, \ \cdots, \ \eta_{\hat{H}_2}, \ \xi_1, \ \cdots, \ \xi_{\hat H_1})$
be
\begin{equation}
\label{eqn:min2}
\sum\limits_{i=1}^m \sum_{t=1}^T
\Big(
y_{it} - {\bf B}^{\rm T}({\bf X}_{it}^{\rm T}\mbox{\boldmath $\beta$}_i)\mbox{\boldmath $\theta$}_i
\Big)^2.
\end{equation}
with $\beta_{ij}$, $i=1, \ \cdots, \ m$, $j=1, \ \cdots, \ p$, being replaced
by $\eta_k$ if $(i, j) \in \hat{G}_{2, k}$, and $\theta_{ij}$, $i=1, \ \cdots,
\ m$, $j=1, \ \cdots, \ K$, being replaced by $\xi_s$ if
$(i, j) \in \hat G_{1,s}$.  Let
$(\hat{\eta}_1, \ \cdots, \ \hat{\eta}_{\hat{H}_2}, \ \hat{\xi}_1, \ \cdots, \
\hat{\xi}_{\hat H_1})$ minimise
$L(\eta_1, \ \cdots, \ \eta_{\hat{H}_2}, \ \xi_1, \ \cdots, \ \xi_{\hat H_1})$.
The final estimator $\hat{\beta}_{ij}$ of $\beta_{ij}$ is $\hat{\eta}_k$ if
$(i, j) \in \hat{G}_{2, k}$, and the final estimator $\hat{\theta}_{ij}$ of
$\theta_{ij}$ is $\hat{\xi}_s$ if $(i, j) \in \hat G_{1,s}$.  Once we have the
estimator $\hat{\theta}_{ij}$, the estimator $\hat{g}_i(\cdot)$ of $g_i(\cdot)$
is taken to be ${\bf B}(\cdot)^{{\mbox{\rm\tiny T}}} \hat{\mbox{\boldmath $\theta$}}_i$.

\end{enumerate}

\bigskip

\noindent
\begin{rmk}  When dealing with the unknown functions $g_i(\cdot)$,
$i=1, \ \cdots, \ m$, in the estimation procedure, instead of treating each
unknown function as a single undivided unit to conduct homogeneity pursuit, we
work on the coefficients of its B-Spline decomposition.  This is because there
may still be some kind of homogeneity between two functions even if they are
different.  For example, for two different functions, it could be the case that
some coefficients of the B-Spline decomposition of one function are the same as
some coefficients of the B-Spline decomposition of another one.  If we treat
each unknown function as a single undivided unit to conduct homogeneity
pursuit, we would not identify or use this kind of homogeneity, which would
make our final estimators not as efficient as they should.
\end{rmk}

\subsection{Computational algorithm}
\label{est2}

In the estimation procedure described in Section \ref{est1}, the minimiser of
(\ref{eqn:min1}) does not have a closed form, neither does the minimiser of
$L(\eta_1, \ \cdots, \ \eta_{\hat{H}_2}, \ \xi_1, \ \cdots, \ \xi_{\hat H_1})$.  To
conduct the minimisation of either of the two objective functions, we appeal to
the standard NLS algorithm, and use the {\tt nlsLM of minpack.lm} package in R
to implement it.   One
can also use other NLS software, for example, the NLS routine {\tt lsqnonlin()}
from MATLAB and {\tt PROC NLIN} from SAS.  To use the {\tt nlsLM of minpack.lm}
package in R, we first need to find an initial value.  The initial value for
minimising (\ref{eqn:min1}) can be obtained as follows:
\begin{enumerate}
\item[(1)] Apply the standard least squares estimation for the linear models to
$(y_{it}, \ {\bf X}_{it})$, $t = 1, \ \cdots, \ T$, and denote the resulting
estimator by $\check{\mbox{\boldmath $\beta$}}_i$, the initial value for $\mbox{\boldmath $\beta$}_i$ is taken to
be $\mbox{\boldmath $\beta$}_i^{(0)} = \check{\beta}_{i0}^{-1} \check{\mbox{\boldmath $\beta$}}_i$,
$\check{\beta}_{i0}$ is the first component of $\check{\mbox{\boldmath $\beta$}}_i$.

\item[(2)] Substitute $\mbox{\boldmath $\beta$}_i^{(0)}$ for $\mbox{\boldmath $\beta$}_i$ in (\ref{eqn:min1}),
then minimise (\ref{eqn:min1}) with respect to $\mbox{\boldmath $\theta$}_i$, the minimiser
$\mbox{\boldmath $\theta$}_i^{(0)}$ is the initial value of $\mbox{\boldmath $\theta$}_i$.
\end{enumerate}
Once we have $\mbox{\boldmath $\beta$}_i^{(0)}$ and $\mbox{\boldmath $\theta$}_i^{(0)}$, the minimiser of
(\ref{eqn:min1}) can be obtained by the {\tt nlsLM of minpack.lm} package in R
straightforwardly.

For any set $A$, let $|A|$ be the number of elements in $A$.  The initial
value for minimising $L(\eta_1, \ \cdots, \ \eta_{\hat{H}_2}, \ \xi_1, \
\cdots, \ \xi_{\hat H_1})$ can
be obtained through the initial estimates of $\mbox{\boldmath $\beta$}_i$ and $\mbox{\boldmath $\theta$}_i$,
obtained in Stage 1 of the estimation procedure in Section \ref{est1}, as
follows:
$$
\eta_s^{(0)}
=
\left(|\hat{G}_{2,s}| \right)^{-1}
\sum\limits_{(i,j) \in \hat{G}_{2,s}} \tilde{\beta}_{ij},
\quad
s = 1, \ \cdots, \ \hat{H}_2
$$
and
$$
\xi_s^{(0)}
=
\left(|\hat G_{1,s}| \right)^{-1}
\sum\limits_{(i,j) \in \hat G_{1,s}} \tilde{\theta}_{ij},
\quad
s = 1, \ \cdots, \ \hat H_1.
$$
Once we have the initial value $(\eta_1^{(0)}, \ \cdots, \
\eta_{\hat{H}_2}^{(0)}, \ \xi_1^{(0)}, \ \cdots, \ \xi_{\hat H_1}^{(0)})$, we
can have the minimiser of $L(\eta_1, \ \cdots, \ \eta_{\hat{H}_2}, \
\xi_1, \ \cdots, \ \xi_{\hat H_1})$ by using the {\tt nlsLM of minpack.lm}
package in R straightforwardly.

\subsection{Selection of tuning parameters}
\label{tune}

The threshold $\delta$ in the Stage 2 of the proposed estimation procedure,
described in Section \ref{est1}, plays a key role for the success of the
homogeneity pursuit.  As far as the implementation of the homogeneity pursuit
is concerned, the selection of $\delta$ is equivalent to the selection of
$\hat H_1$ and $\hat H_2$, and to select an integer is easier, therefore,
in this section, instead of selecting $\delta$, we develop a cross-validation
procedure to select the two tuning parameters, $\hat H_1$ and $\hat H_2$.

For the single index model (\ref{model}) where $\mbox{\bf X}_{it}$'s are independent
across $t=1, \ \cdots, \ T,$  we implement a $L$-fold cross validation
approach.  In particular, for a given pair $\{H_1, \ H_2\}$, we remove
$1/L$th of the observed time points for
$\{(y_{it}, \ {\bf X}_{it}), \ i=1, \ \cdots, \ m, \ t=1, \ \cdots, \ T\}$
as a validation set, estimate the single index model (\ref{model}) with
identified homogeneity structure on the remaining data, compute the squared
error between $y_{it}$ and fitted values
$\hat g_i({\bf X}_{it}^{{\mbox{\rm\tiny T}}}\hat\mbox{\boldmath $\beta$}_i)
={\bf B}^{\rm T}({\bf X}_{it}^{\rm T}\hat\mbox{\boldmath $\beta$}_i)\hat\mbox{\boldmath $\theta$}_i$, on the validation set,
and repeat this procedure $L$ times to calculate the cross-validated mean
squared error and its corresponding standard error.  We search over a grid of
$\{H_1, \ H_2\}$ values and apply the one-standard-error rule to choose the
smallest model for which the estimated cross-validated error is within one
standard error of the lowest point on the error surface.  The rationale here is
that if a set of models appear to be more or less equally good, then we might
tend to choose the simplest model.  Across the candidate pairs,
$\{\hat H_1, \ \hat H_2\},$ whose corresponding errors are within this
deviation, one can choose the smallest $\hat H_1$ after selecting the smallest
$\hat H_2$ or switch the selection order or select the smallest value of
$\hat H_1+\hat H_2$, we take the first approach since it produces better model
selection consistency in our numerical experiments.  A similar
one-standard-deviation-rule technique has been adopted to choose the
regularisation parameter with a smaller model size for the lasso problems
\cite[]{james13}.

When $\mbox{\bf X}_{it}$'s are time dependent panel data, we implement a rolling
procedure to perform cross-validation for time series.  More specifically,
for each $r=L, \ L-1, \ \cdots, \ 1,$ we rollingly treat
$\{(y_{it}, \ {\bf X}_{it}), \ i=1, \ \cdots, \ m, \ t=1, \ \cdots, \ T-r\}$ as
training observations and $\{(y_{i,T-r+1}, \ {\bf X}_{i,T-r+1}), \ i=1, \ \cdots, \
m\}$ as validation set, calculate the squared error between each $y_{it}$ and
its fitted value.  Finally, we apply the one-standard-deviation-rule on the
lowest cross-validated mean squared error and choose $\hat H_2$ and $\hat H_1.$

In the cross-validation procedure when we need make predictions for validation
set, the domain in ${\bf B}({\bf X}_{it}^{{\mbox{\rm\tiny T}}}\hat \mbox{\boldmath $\beta$}_i)$ for traning data set might
not cover that for validation set.  We adopt the idea in \cite{wangyang09} by
mapping ${\bf X}^{{\mbox{\rm\tiny T}}}_{it}\mbox{\boldmath $\beta$}_i$ to
$F_i({\bf X}^{{\mbox{\rm\tiny T}}}_{it}\mbox{\boldmath $\beta$}_i) \sim \text{Unif}[0, \ 1]$,
where $F_i$ is the distribution function of ${\bf X}^{{\mbox{\rm\tiny T}}}_{it}\mbox{\boldmath $\beta$}_i$.  We then
implement the estimation procedure described in Section~\ref{est1} by
decomposing
$g_i({\bf X}^{{\mbox{\rm\tiny T}}}_{it}\mbox{\boldmath $\beta$}_i)
\approx
{\bf B}^{{\mbox{\rm\tiny T}}}\big(F_i({\bf X}^{{\mbox{\rm\tiny T}}}_{it}\mbox{\boldmath $\beta$}_i)\big) \mbox{\boldmath $\theta$}_i$.
The proposed approach is thus able to make predictions and, as demonstrated by
some numerical studies, provides very similar sample performance in terms of
estimation accuracy.

\subsection{Post-processing step}
\label{post}
We equip the Binary-Segmentation-algorithm-based homogeneity pursuit with an
additional step aimed to enhance the accuracy of detected change-points
locations through a fine-scale search. To be specific, at each change-point,
we re-calculate  $\Delta_{ij}(\kappa)$ over the interval between two adjacent
change-points and identify the new change-point location to replace the old
one. We perform this post-processing procedure by iteratively cycling
through all neighbouring change-points and fine-tuning the change-points
locations. This procedure is terminated when the set of change-points does not
change. Our numerical experiments show that this extra post-processing step
apparently improve the accuracy of each estimated change-point location and
hence the identified homogeneity structure for model~\eqref{model}.

\section{Asymptotic properties}
\label{asy}

In this section, we are going to investigate the asymptotic behaviour of the
estimators obtained by the proposed estimation procedure, which we call
correct-fitting, and compare with the estimators obtained without homogeneity
pursuit, which is the initial estimators obtained in the Stage 1 in the
proposed estimation procedure, we call it over-fitting, and the estimators
obtained under the assumption that all individuals share the same index
(namely, $\mbox{\boldmath $\beta$}_1 = \cdots = \mbox{\boldmath $\beta$}_m$), which we call under-fitting.  The
asymptotic theory presented in this section is in the sense that
$T\longrightarrow \infty$, and $m$, $p$ are all possibly diverging to infinity
but $H_1$, $H_2$ are fixed.  This agrees with many applications in which $H_1$
and $H_2$ are expected to be small and thus significant reduction of unknown
parameters can be achieved by clustering the parameters.  To make the
presentation neat, we state the asymptotic theorems in this section and leave
all technical proofs in the Appendix.

Let $\epsilon_{t}=(\epsilon_{1t}, \ \cdots, \ \epsilon_{mt})^{\rm T}$,
$$
\overline{\bf X}_{it}=(X_{it,1}, \ \cdots, \ X_{it,p})^{\rm T},
\quad
y_t=(y_{1t}, \ \cdots, \ y_{mt})^{\rm T},
\quad
{\bf X}_{t}=({\bf X}_{1t}^{\rm T}, \ \cdots, \ {\bf X}_{mt}^{\rm T})^{\rm T}.
$$
In this paper, we assume $(y_{t}, \ {\bf X}_{t}, \ \epsilon_t)$ are stationary
with $\alpha(l)\le \rho^l$ for some $\rho<1$, and $\epsilon_t$ is
independent of ${\bf X}_t$.  Note that unlike \cite{vogt15}, we do not need to
assume independence or stationarity of variables cross $i$.

We start with the asymptotic properties of the estimators obtained without
homogeneity pursuit. The convergence rate of the estimator $\tilde{\beta}_{ij}$ is of order $T^{-1/2}$, and the convergence rate of the estimator
$\tilde{g}_i(u)$ is of order $T^{-2/5}$, which is as expected as we assumed the
functions are twice differentiable.

\bigskip

\begin{thm}\label{thm:new1}
({\it Over-fitting case}).  For any $i$, $i=1, \ \cdots, \ m$, and
$1\le j\le p$, under the conditions (C1)-(C4) and (C5') in the Appendix, we
have
$$
T^{1/2} (\widetilde\varepsilon_{ij}^{\rm T}\widetilde\mbox{\boldmath $\Theta$}_2\widetilde\varepsilon_{ij})^{-1/2}
\left(\tilde{\beta}_{ij} - \beta_{ij} \right)
\stackrel{D}{\longrightarrow} N(0,\ 1)
$$
and
$$
T^{2/5}(\widetilde{\bf b}_i^{\rm T}(u)\widetilde\mbox{\boldmath $\Theta$}_1\widetilde{\bf b}_i(u))^{-1/2}
\left(\widetilde{g}_i(u) - g_{i}(u)-r_i(u) \right)
\stackrel{D}{\longrightarrow} N(0,\ 1),
$$
where $\widetilde\varepsilon_{ij}$ and $\widetilde{\bf b}_i(u)$ are unit vectors,  $\widetilde\mbox{\boldmath $\Theta$}_1$,
$\widetilde\mbox{\boldmath $\Theta$}_2$ are matrices with eigenvalues bounded and bounded away from
zero, all these quantities are defined in the proof in the
Appendix \ref{sec:conv4}.  The bias term
$r_i(u)=g_{i}(u)-{\bf B}^{\rm T}(u)\mbox{\boldmath $\theta$}_{0i}$ satisfies $|r_i(u)|\le CK^{-2}$,
where $\mbox{\boldmath $\theta$}_{0i}$ is the vector of spline coefficients used to approximate
$g_i$ as defined in Appendix \ref{sec:ass}
parameter of $g_i$ as defined in assumption  (C3).
\end{thm}

\bigskip

Let $m_i$ be the size of $G_{1,h}$ that contains $i$, and $m_{ij}$ be the size
of $G_{2,h}$ that contains $\beta_{ij}$.  To make the statement about the
correct-fitting case cleaner, we assume that all $m_i$ are of the same order
and all $m_{ij}$ are of the same order ($\max_{i,j} m_{ij}/\min_{i,j} m_{ij}$
and $\max_i m_i/\min_i m_i$ are bounded) in the following theorem, which shows
in particular that the convergence rate of the estimator $\hat{\beta}_{ij}$ is
of order $(mpT)^{-1/2}$, and the convergence rate of the estimator
$\hat{g}_i(u)$ is of order $(mT)^{-2/5}$.

\bigskip

\begin{thm}\label{thm:new2} ({\it Correct-fitting case}).  For any $i$,
$i=1, \ \cdots, \ m$, and $1\le j\le p$, under the conditions (C1)-(C6) in the
Appendix, we have
$$
(mpT)^{1/2} (\varepsilon_{ij}^{\rm T}\mbox{\boldmath $\Theta$}_2\varepsilon_{ij})^{-1/2}
\left(\hat{\beta}_{ij} - \beta_{ij} \right)
\stackrel{D}{\longrightarrow} N(0, \ 1)
$$
and
$$
(mT)^{2/5}({\bf b}_i^{\rm T}(u)\mbox{\boldmath $\Theta$}_1{\bf b}_i(u))^{-1/2}
\left(\hat{g}_i(u) - g_i(u)-r_i(u) \right)
\stackrel{D}{\longrightarrow} N(0, \ 1),
$$
where $\varepsilon_{ij}$ and ${\bf b}_i(u)$ are unit vectors, $\mbox{\boldmath $\Theta$}_1$, $\mbox{\boldmath $\Theta$}_2$ are
matrices with eigenvalues bounded and bounded away from zero, all these
quantities are defined in the proof in the Appendix \ref{sec:conv4}.
\end{thm}

\bigskip

Finally, for the under-fitting case, let $\check{\mbox{\boldmath $\beta$}}_i$ and
$\check{g}_i(\cdot)$ be the estimators of $\mbox{\boldmath $\beta$}_i$ and $g_i(\cdot)$ obtained
under the assumption that all individuals share the same unknown parameters.

\bigskip

\begin{thm}\label{thm:new3} ({\it Under-fitting case}).  Suppose the
$\mbox{\boldmath $\beta$}_i$s are sufficiently separated in the sense that for
$\bar\mbox{\boldmath $\beta$}:=\sum_{i=1}^m\mbox{\boldmath $\beta$}_i/m$,
$$\frac{1}{mp}\sum_{i=1}^m\|\mbox{\boldmath $\beta$}_i-\bar\mbox{\boldmath $\beta$}\|^2\ge c$$
for some $c>0$, then
$$\frac{1}{mp}\sum_{i=1}^m\|\check \mbox{\boldmath $\beta$}_i-\mbox{\boldmath $\beta$}_i\|^2\ge c.$$
Similarly, if
$$\frac{1}{m}\sum_{i=1}^m\int\left |{g}_i(u) - \bar g(u) \right|^2 du \ge c,
$$
where $\bar g(u)=m^{-1}\sum_{i=1}^m g_i(u)$, then
$$\frac{1}{m}\sum_{i=1}^m\int\left |\check{g}_i(u) -  g_i(u) \right|^2 du \ge c.
$$
\end{thm}

\section{Simulation studies}
\label{simu}

In this section, we are going to use a simulated example to demonstrate how
accurate the proposed estimation is.  We will also show much loss it would
inflict if the homogeneity structure is ignored or mistakenly specified as that
all individuals share the same index coefficients or the same link function.

\noindent
{\bf Example}.  We generate a sample from model (\ref{model}) with $p=2$ and
WLOG an even $m,$
$$
g_i(u)
=
\begin{cases}
\sin(\pi u/4) & \text{when } i=1, \ 2, \ \cdots, \ m/2,
\\
\cos(\pi u/4) & \text{when } i=m/2+1, \ \cdots, \ m,
\end{cases}
$$
and
$$
\mbox{\boldmath $\beta$}_i
=
\begin{cases}
(1, \ -1.5\sqrt{0.2}, \ -0.5\sqrt{0.2})^{{\mbox{\rm\tiny T}}} &
\text{when } i=1, \ 3, \ \cdots, \ m-1,
\\
(1, \ 0.5\sqrt{0.2}, \ 1.5\sqrt{0.2})^{{\mbox{\rm\tiny T}}} &
\text{when } i=2, \ 4, \ \cdots, \ m,
\end{cases}
$$
where $\|\mbox{\boldmath $\beta$}_i\|^2=1.5$ for $i=1, \ \cdots, \ m$.  Let ${\bf X}_{it}$ and
$\epsilon_{it}$, $i=1, \ \cdots, \ m$, $t=1, \ \cdots, \ T$ be independently
generated from $\frac{1}{\sqrt{1.5}}N(\mbox{\bf 0}_3, \ {\bf I}_3)$ truncated by
$[-1.343, \ 1.343]^3$ (the range of 5th to 95th quantiles for $N(0, \ 2/3)$)
and $N(0, \ \sigma^2)$, respectively.  Once ${\bf X}_{it}$ and $\epsilon_{it}$ are
generated, $y_{it}$ can be generated through (\ref{model}).

We conduct the simulated example for various $m$s and $T$s with $\sigma=0.2$,
and compare our proposed approach to its potential competitors based on the
following performance metrics:
\begin{enumerate}
\item[(1)]{\it Estimation accuracy}.  For an estimator $\hat{\mbox{\boldmath $\beta$}}_i$ of
$\mbox{\boldmath $\beta$}_i$, we use the mean squared error (MSE), namely
$\mbox{MSE}(\hat{\mbox{\boldmath $\beta$}}_i) = E\left(\|\hat{\mbox{\boldmath $\beta$}}_i - \mbox{\boldmath $\beta$}_i \|^2\right)$, to
assess the estimation
error of $\hat{\mbox{\boldmath $\beta$}}_i$.  Analogously, for an estimator $\hat g_i(\cdot)$ of
$g_i(\cdot)$, its estimation accuracy can be evaluated based on the mean
integrated squared error,
$$
\text{MISE}(\hat g_i)
=
E \left\{\int\left(\hat{g}_i(u) - g_i(u)\right)^2 du \right\}.
$$
To avoid the situation where the performance is dominated by the poor boundary
behaviour, we let the integral domain to be non-boundary region, which is
between the 1st and 99th quantiles of
$\{{\bf X}_{it}^{{\mbox{\rm\tiny T}}}\mbox{\boldmath $\beta$}_i, \ t=1, \ \cdots, \ T\}$.

\item[(2)]{\it Homogeneity structure identification consistency}.  To evaluate
the distance between the detected homogeneity structure and the true one, we
use the normalized mutual information (NMI) \cite[]{ke15}, which measures the
similarity between two partitions.  Suppose
$\mathbb{C}=\{C_1, \ C_2, \ \cdots \}$ and
$\mathbb{D}=\{D_1, \ D_2, \ \cdots \}$ are two partitions of
$\{1, \ \cdots, \ n\}$, the NMI is defined as
$$
\text{NMI}(\mathbb{C}, \ \mathbb{D})
=
\frac{I(\mathbb{C}, \ \mathbb{D})}
{[H(\mathbb{C})+ H(\mathbb{D})]/2},
$$
where
$$
I(\mathbb{C}, \ \mathbb{D})
=
\sum\limits_{k,j}
\big(|C_k \cap D_j|/n\big) \hbox{log}\big(n| C_k \cap D_j|/|C_k||D_j|\big)
$$
and
$$
H(\mathbb{C})
=
-\sum\limits_{k}
\big(|C_k|/n\big) \hbox{log}\big(|C_k|/n\big).
$$
The NMI takes values in $[0, \ 1]$ with larger values indicating higher level
of similarity between two partitions.  For an estimated partition
$\hat{\mathbb{G}}_2=\{\hat G_{2,1}, \ \cdots, \ \hat G_{2, \hat H_2}\}$ of
$\{(i, \ j): \ 1, \ \cdots, \ m, \ j=1, \ \cdots, \ p\}$, obtained in the
Stage 2 of the proposed estimation procedure in Section \ref{est1}, we
calculate $\text{NMI}(\hat{\mathbb{G}}_2, \ \mathbb G_2)$ to assess how close
to the true homogeneity structure in $\beta_{ij}$s the estimated one is.
Similarly, for an estimated partition $\hat{\mathbb{G}}_1$ of
$\{i: \ 1, \ \cdots, \ m\}$, we use
$\text{NMI}(\hat{\mathbb{G}}_1, \ \mathbb G_1)$ to evaluate how close
the estimated homogeneity structure in $g_{i}(\cdot)$s is to the true one.
\end{enumerate}

For each case, we apply either the single index model (\ref{model}) with the
standard estimation procedure, the initial estimation of the proposed
estimation procedure in Section~\ref{est1}, which we call over-fitting (Over),
the single index model (\ref{model}) with the homogeneity structure
(\ref{con}) together with the proposed estimation procedure, which we call
correct-fitting, the single index model (\ref{model}) with all individuals
share the same index vector (namely, $\mbox{\boldmath $\beta$}_1 = \cdots = \mbox{\boldmath $\beta$}_m$),
which we call Under-I, the single index model (\ref{model}) with all
individuals share the same link function (namely,
$g_1(\cdot) = \cdots = g_m(\cdot)$, i.e. $\mbox{\boldmath $\theta$}_1 = \cdots = \mbox{\boldmath $\theta$}_m$),
which we call Under-F, or the single index model (\ref{model}) with all
individuals share both the same index and link function, which we call
Under-I-F, to the simulated data set.

We develop three methods under the correct-fitting case.  The first approach,
named Correct-C, optimises (\ref{eqn:min2}) based on the estimated
{\it componentwise} homogeneity structure in $\beta_{ij}$s and $\theta_{ij}$s,
obtained in the Stage 2 of the proposed estimation procedure in Section
\ref{est1} with the tuning parameters selected through the cross-validation
approach described in Section~\ref{tune}.  The second approach, Correct-V, is
the same as the first approach but optimises (\ref{eqn:min2}) based on the
estimated componentwise homogeneity structure in $\beta_{ij}$s and
{\it vectorwise} homogeneity structure in $\mbox{\boldmath $\theta$}_i$s which can be obtained
through the estimated componentwise homogeneity structure in $\theta_{ij}$s.
The third approach, which we call Correct-NMI, is the same as the second
approach but with the tuning parameters selected to be the one maximising
$\text{NMI}(\hat{\mathbb{G}}_2, \ \mathbb G_2)$ and
$\text{NMI}(\hat{\mathbb{G}}_1, \ \mathbb G_1)$.  In practice without knowing
the true homogeneity structure, one cannot implement Correct-NMI.

Under-I, Under-F and Under-I-F are three kinds of under-fitting.  For Under-I
or Under-F, the homogeneity structure in $\mbox{\boldmath $\theta$}_i$s or $\beta_{ij}$s is
estimated in the same way as that in the proposed estimation
procedure in Section \ref{est1} with the tuning parameters still selected by
the one-standard-deviation-rule cross-validation approach.

We compare over-fitting, correct-fittings and under-fittings to the oracle case
where the true homogeneity structure is used.  The computational algorithms
for the under-fitting and oracle estimators are the same as that for the
correct-fitting, but use either identified or pre-specified homogeneity
structure.  We compare the sample performance of all eight approaches in our
conducted simulation study.

We report the results for estimation errors and NMIs for $\mbox{\boldmath $\beta$}_i$s and
$g_i(\cdot)$s averaged over 100 replicates in Tables~\ref{tab_beta} and
\ref{tab_fun}, respectively.  In terms of estimation error, the overall
estimation accuracy is improved as $m$ and $T$ increase and three
correct-fitting approaches perform very well as reflected in their lower values
of MSEs and MISEs.  Among the three methods, Correct-NMI provides the best
performance even producing very comparable MSEs and MISEs with the oracle
estimator and Correct-C is outperformed by Correct-V in most settings.  This is
somewhat expected, since, unlike Correct-C, which optimises (\ref{eqn:min2})
based on the detected homogeneity structure in $\beta_{ij}$s and
$\theta_{ij}$s, Correct-V separates the final estimation step from the
cross-validation procedure, which is used to identify the homogeneity structure
in $\beta_{ij}$s and $\mbox{\boldmath $\theta$}_{i}$s.   Analogously, Correct-NMI solves a
separate optimisation after detecting the homogeneity structure based on the
largest NMIs.  It is also worth noting that the over-fitting and under-fitting
methods, which either ignores or mistakenly specify the homogeneity structure,
provide much worse results, highlighting the importance of incorporating the
appropriate homogeneity structure.  In terms of selecting the structure of
homogeneity, we observe that three correct-fitting methods produce perfect
identifications of homogeneity structure in $\beta_{ij}$s and Correct-NMI
provides the largest NMI values indicating that it can effectively recover the
true homogeneity structure in $g_i(\cdot)$s.  The performance of Correct-C and
Correct-V deteriorates when $m$ increases, this is intuitively due to the
increased $m$ values and the cross-validation procedure, which tends to choose
a larger number of change points as $m$ increases, resulting in smaller NMI
values for $\hat{\mathbb{G}}_1$.

\begin{table}[htbp]
	\begin{center}
		\caption{{\bf The Average of $\mbox{MSE}(\hat{\mbox{\boldmath $\beta$}}_i)$,
$i=1, \ \cdots, \ m$, and Average NMIs for $\hat{\mathbb{G}}_2$}
\centerline{{\it All entries for MSEs are $10^4$ times their actual values}}}
		\label{tab_beta}
		\vspace{0.3cm}
		\begin{tabular}{cc|ccc|ccc}
			\hline
			\hline
			& $T$ & \multicolumn{3}{c|}{400} & \multicolumn{3}{c}{800}\tabularnewline
			& $m$ & 30 & 60 & 90 & 30 & 60 & 90\tabularnewline
			\hline
			\multirow{8}{*}{MSE}
			& Oracle & 0.380 & 0.205 & 0.133 & 0.214 & 0.128 & 0.101\tabularnewline
			& Correct-C & 0.473 & 0.221 & 0.147 & 0.157 & 0.093 & 0.085 \tabularnewline
			& Correct-V & 0.381 & 0.208 & 0.134 & 0.214 & 0.129 & 0.100 \tabularnewline
			& Correct-NMI & 0.381 & 0.206 & 0.133 & 0.214 & 0.128 & 0.101 \tabularnewline
			& Over & 5.438 & 5.285 & 5.246 & 2.636 & 2.689 & 2.663 \tabularnewline
			& Under-I-F & 4005.6 & 4002.6 & 4002.1 & 4003.4 & 4001.5 & 4001.1\tabularnewline
			& Under-I & 4002.4 & 4001.7 & 4001.3 & 4001.8 & 4001.0 & 4000.6\tabularnewline
			& Under-F & 21.233 & 3.393 & 2.259 & 3.111 & 1.454 & 1.123\tabularnewline
			\hline
			\multirow{8}{*}{NMI}
			& Oracle & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000\tabularnewline
			& Correct-C & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000\tabularnewline
			& Correct-V & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000\tabularnewline
			& Correct-NMI & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000\tabularnewline
			& Over & 0.339 & 0.290 & 0.267 & 0.339 & 0.290 & 0.267\tabularnewline
			& Under-I-F & 0 & 0 & 0 & 0 & 0 & 0\tabularnewline
			& Under-I & 0 & 0 & 0 & 0 & 0 & 0\tabularnewline
			& Under-F & 0.817  & 0.810 & 0.807 & 0.812 & 0.809 & 0.806\tabularnewline
			\hline
			\hline
		\end{tabular}
	\end{center}
\end{table}

\begin{table}[htbp]
	\begin{center}
		\caption{{\bf The Average of $\mbox{MISE}(\hat{g}_i)$, $i=1, \
\cdots, \ m$, and Average NMIs for $\hat{\mathbb{G}}_1$}
\centerline{{\it
where all entries for MISEs are $10^2$ times their actual values
}}}
		\label{tab_fun}
				\vspace{0.3cm}
		\begin{tabular}{cc|ccc|ccc}
			\hline
			\hline
			& $T$ & \multicolumn{3}{c|}{400} & \multicolumn{3}{c}{800}\tabularnewline
			& $m$ & 30 & 60 & 90 & 30 & 60 & 90\tabularnewline
			\hline
			\multirow{8}{*}{MISE}
			& Oracle & 0.260 & 0.249 & 0.243 & 0.251 & 0.241 & 0.236\tabularnewline
			& Correct-C & 0.667 & 0.553 & 0.485 & 0.366 & 0.318 & 0.307 \tabularnewline
			& Correct-V & 0.300 & 0.289 & 0.311 & 0.257 & 0.246 & 0.246\tabularnewline
			& Correct-NMI & 0.266 & 0.254 & 0.249 & 0.252 & 0.241 & 0.238\tabularnewline
			& Over & 0.548 & 0.544 & 0.540 & 0.410 & 0.403 & 0.396\tabularnewline
			& Under-I-F & 87.988 & 87.963 & 87.959 & 88.746 & 88.726 & 88.647\tabularnewline
			& Under-I & 10.099 & 9.853 & 9.952 & 9.770 & 9.719 & 9.711\tabularnewline
			& Under-F & 85.196 & 85.240 & 85.235 & 86.017 & 86.035 & 85.952\tabularnewline
			\hline
			\multirow{8}{*}{NMI}
			& Oracle & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \tabularnewline
			& Correct-C & 0.816 & 0.790 & 0.637 & 0.958 & 0.948 &0.881 \tabularnewline
			& Correct-V & 0.816 & 0.790 & 0.637 & 0.958 & 0.948 & 0.881\tabularnewline
			& Correct-NMI & 0.965 & 0.966 & 0.948 & 0.994 & 0.997 & 0.997\tabularnewline
			& Over & 0.339 & 0.290 & 0.267 & 0.339 & 0.290 & 0.267\tabularnewline
			& Under-I-F & 0 & 0 & 0 & 0 & 0 & 0\tabularnewline
			& Under-I & 0.811 & 0.797 & 0.783 & 0.901 & 0.894 & 0.900\tabularnewline
			& Under-F & 0 & 0 & 0 & 0 & 0 & 0\tabularnewline
			\hline
			\hline
		\end{tabular}
	\end{center}
\end{table}

\section{Real data analysis}
\label{re}

We will illustrate the proposed method with two real data examples in this
section.

\subsection{Industrial Portfolio's return}
\label{port}

We first study the data set about $m=49$ Industrial Portfolios' daily simple
return from 1/8/2015 to 31/12/2015.  This data set can be freely downloaded
from Kenneth French's website

\begin{center}
{\tt http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data\_library.html}
\end{center}

This data set has been analysed in quite a few literature.  For example,
\cite{guo2016dynamic} used this data set to demonstrate the performance of a
newly developed dynamic portfolio allocation.  In this paper, we are going to
explore the homogeneity structure in this data set by our proposed method.

Let $y_{it}$ be the daily simple return of the $i$th portfolio at the $t$th
day, $i=1, \ \cdots, \ m, \ t=1, \ \cdots, \ T$, and
${\bf X}_{it}=(X_{t1}, \ X_{t2}, \ X_{t3})^{{\mbox{\rm\tiny T}}}$ be the observation of the
Fama-French three factors, where $X_{t1}, \ X_{t2}, \ X_{t3}$ respectively
represent the market (Rm-Rf), size (SMB) and value (HML) factors at the $t$th
day.

We apply the single index model (\ref{model}) with the unknown homogeneity
structure (\ref{con}) to fit the data set.  From interpretation point of view,
the homogeneity structure in the unknown link functions, $g_i(\cdot)$s, where
each $g_i(\cdot)$ is treated as a single undivided unit would make much more
sense than the homogeneity structure in the coefficients of the B-Spline
decompositions of $g_i(\cdot)$s.  Therefore, we use the Correct-V, described
in Section~\ref{simu}, to identify the homogeneity structure in $\beta_{ij}$s
or $g_i(\cdot)$s, and estimate the unknown parameters and unknown functions.

In the implementation of the Correct-V, we implement the method in
Section \ref{est1} with the tuning parameters selected by the cross-validation
for time series as described in Section~\ref{tune}.  Specifically, we define
the cross-validated mean squared error
\begin{equation}
\label{cv.err}
\mbox{CV}
=
\frac{1}{m L}\sum\limits_{i=1}^{m} \sum\limits_{t=T-L+1}^{T}
\left(y_{it}-\hat{y}_{it}\right)^2
\end{equation}
where $L=30$.  Note that we here do not apply the one-standard-rule when
performing the cross
validation to select the tuning parameters for identifying the homogeneity
structure, since we have already selected a small enough model with 11 and 2
detected groups in index coefficients and link functions, respectively.
Table~\ref{tab4} provides the identified clustering results for $\beta_{i2}$,
$\beta_{i3}$, $g_{i}(\cdot)$, $i=1, \ \cdots, \ 49$, and
Figure~\ref{funcplot_port} plots the estimated link functions.  We observe a
few apparent patterns.  Firstly, the estimated link functions are very linear
indicating the linear relationship between portfolio returns and Fama-French
three factors, which has been verified by broad empirical studies.  Secondly,
many portfolios belonging to similar industrials were grouped into the same
cluster for the estimated index coefficients, e.g. Hardw, Softw and Agric,
Food, Soda were clustered into Groups~8 and 4 in terms of the estimated
coefficients for factors SMB and HML respectively.

\begin{table}[htbp]
\footnotesize
	\caption{\label{tab4}{\bf
Grouping Results for The Index Coefficients for SMB, HML and
\centerline{Link Functions of 49 Industrial Portfolios}}}
\begin{center}
\begin{tabular}{c|cccccccccc}
	\hline
	\hline
	& Agric & Food & Soda & Beer & Smoke & Toys & Fun & Books & Hshld & Clths\tabularnewline
	SMB & 6 & 5 & 6 & 4 & 2 & 8 & 9 & 7 & 6 & 6\tabularnewline
	HML & 4 & 4 & 4 & 3 & 4 & 5 & 3 & 4 & 5 & 6\tabularnewline
	Function & ii & i & ii & ii & i & i & i & i & i & ii\tabularnewline
	\hline
	& Hth & MedEq & Drugs & Chems & Rubbr & Txtls & BldMt & Cnstr & Steel & FabPr\tabularnewline
	SMB & 9 & 8 & 9 & 7 & 7 & 8 & 8 & 8 & 9 & 10\tabularnewline
	HML & 2 & 2 & 1 & 6 & 4 & 4 & 6 & 6 & 9 & 7\tabularnewline
	Function & i & i & i & i & i & i & i & i & i & ii\tabularnewline
	\hline
	& Mach & ElcEq & Autos & Aero & Ships & Guns & Gold & Mines & Coal & Oil\tabularnewline
	SMB & 8 & 9 & 7 & 6 & 8 & 6 & 10 & 9 & 11 & 9 \tabularnewline
	HML & 8 & 6 & 6 & 6 & 7 & 4 & 10 & 7 & 11 & 9 \tabularnewline
	Function & i & i & i & i & i & i & ii & i & ii & i\tabularnewline
	\hline
	& Util & Telcm & PerSv & BusSv & Hardw & Softw & Chips & LabEq & Paper & Boxes\tabularnewline
	SMB & 4 & 8 & 8 & 7 & 8 & 8 & 7 & 7 & 6 & 6 \tabularnewline
	HML & 5 & 4 & 5 & 4 & 5 & 3 & 5 & 5 & 5 & 6\tabularnewline
	Function & i & i & i & i & i & i & i & i & i & i \tabularnewline
	\hline
	& Trans & Whlsl & Rtail & Meals & Banks & Insur & RIEst & Fin & Other & \tabularnewline
	SMB & 7 & 7 & 7 & 7 & 7 & 6 & 7 & 6 & 7 & \tabularnewline
	HML & 6 & 6 & 4 & 3 & 6 & 5 & 5 & 5 & 6 & \tabularnewline
	Function & i & i & i & i & i & i & i & i & ii & \tabularnewline
	\hline
	\hline
\end{tabular}
\end{center}
\end{table}

\begin{figure*}
	\centering
	\includegraphics[width=16cm,height=16cm]{figure/funcplot1_port.pdf}
	\caption{\label{funcplot_port}{Plots of estimated link functions with
respect to $\{{\bf X}_{it}^{{\mbox{\rm\tiny T}}}\hat \mbox{\boldmath $\beta$}_i, t=1,\dots, T\}$.}}
\end{figure*}

\subsection{UK climate data}
\label{temp}

Our second data set, which is available from the UK Met Office website
\begin{center}
{\tt http://www.metoffice.gov.uk/public/weather/climate-historic},
\end{center}
contains monthly data of the mean daily maximum temperature (TMAX), mean daily
minimum temperature (TMIN), days of air frost (AF), total rainfall (RAIN) and
total sunshine duration (SUN) collected from 37 stations across the UK.  We
first remove the missing values and thus select data during the period of
January 1993 to December 2009 from 16 locations.  We then eliminate the
seasonality and trend effects and standardise the data.  Let $y_{it}$ be the
monthly mean temperature, which can be calculated as (TMAX+TMIN)/2, and
${\bf X}_{it}=(X_{it1}, \ X_{it2}, \ X_{it3})^{{\mbox{\rm\tiny T}}}$ be the observations for AF,
RAIN and SUN, from the $i$th station at the $t$th month, $i=1, \ \cdots, \ 16,
\ t=1, \ \cdots, \ T$,

Like the analysis of the Industrial Portfolio's return data set, we apply the
single index model (\ref{model}) with unknown homogeneity structure (\ref{con})
together with the proposed estimation procedure, Correct-V, to the data set.
Table~\ref{tab6} provides the clustering results for the index coefficients and
link functions, where 4 and 2 groups were selected respectively.
Figure~\ref{funcplot_temp} plots the estimated link functions at 16 stations.

It is very interesting to see, from Table~\ref{tab6}, that Oxford, Hurn,
Eastbourne and Bradford share exactly the same model, which implies the impact
of rainfall or total sunshine duration on monthly mean temperature has exactly
the same pattern in these four areas.  The same finding also appears in the
three areas of Waddington, Sheffield and Heathrow, the two areas of
Ross-On-Wye and Eskdalemuir, and the two areas of Paisley and Leuchars.  If we
only focus on the impact of rainfall on monthly mean temperature, the seven
areas of Waddington, Sheffield, Shawbury, Paisley, Leuchars, Lerwick and
Heathrow would have exactly the same pattern.  Similar finding also appears
for the impact of total sunshine duration on monthly mean temperature.

\begin{table}[htbp]
\footnotesize
	\begin{center}
	\caption{{\bf Grouping Results for The Index Coefficients for RAIN,
SUN and
\centerline{Link Functions at 16 Locations}}}
	\label{tab6}
			\vspace{0.3cm}
		\begin{tabular}{c|cccccccc}
			\hline
			\hline
			& Waddington & Sheffield & Shawbury & Ross-On-Wye & Paisley & Oxford & Leuchars & Lerwick\tabularnewline
			RAIN & 4 & 4 & 4 & 3 & 4 & 3 & 4 & 4 \tabularnewline
			SUN & 1 & 1 & 2 & 2 & 3 & 2 & 3 & 4 \tabularnewline
			Function & i & i & i & ii & i & i & i & i \tabularnewline
			\hline
			& Hurn & Heathrow & Eskdalemuir & Eastbourne & Cambridge & Camborne & Bradford & Armagh \tabularnewline
			RAIN & 3 & 4 & 3 & 3 & 3 & 3 & 3 & 4 \tabularnewline
			SUN & 2 & 1 & 2 & 2 & 1 & 3 & 2 & 2 \tabularnewline
			Function & i & i & ii & i & i & i & i & ii \tabularnewline
			\hline
			\hline
		\end{tabular}
	\end{center}
\end{table}

\begin{figure*}
	\label{fig3}
	\centering
	\includegraphics[width=16cm,height=16cm]{figure/funcplot_temp.pdf}
	\caption{\label{funcplot_temp}{Plots of estimated link functions with
respect to $\{{\bf X}_i^{{\mbox{\rm\tiny T}}}\hat \mbox{\boldmath $\beta$}_{it}, t=1, \dots, T\}$.}}
\end{figure*}

\section{Introduction}

\subsection{Preamble}
\label{pre}

Panel data analysis is an important topic in statistics and econometrics.  The
traditional approach for analysing panel data assumes all individuals share the
same unknown parameters, and uses cluster effects to account for the difference
between individuals.  For example, when the linear models are
used, the coefficients of the covariates are assumed to be the same across all
individuals, i.e.
$$
y_{it} = {\bf X}_{it}^{{\mbox{\rm\tiny T}}} \mbox{\boldmath $\beta$} + \epsilon_{it},
\quad
i=1, \ \cdots, \ m;  \  \  t=1, \ \cdots, \ T,
$$
where $y_{it}$ and ${\bf X}_{it}$, a $(p+1)$-dimensional vector, are respectively
the $t$th observations of the response variable and covariate of the $i$th
individual.  $\epsilon_{it}$, $t=1, \ \ldots, \ T$, are correlated for any
given $i$, and the cluster effects are included in $\epsilon_{it}$.  See
\cite{hsiao14} and the reference therein.  Whilst this modelling idea is
useful when the global trend of the impact of a covariate on the response
variable is of our main interest, it does not tell us anything about the
individual attributes which are sometimes very important.

In order to explore the individual attributes, we need to make them more
concrete and distinctive in modelling.  A simple approach to do so would be
using
\begin{equation}
y_{it} = {\bf X}_{it}^{{\mbox{\rm\tiny T}}} \mbox{\boldmath $\beta$}_i + \epsilon_{it},
\quad
i=1, \ \cdots, \ m;  \  \  t=1, \ \cdots, \ T,
\label{eq1}
\end{equation}
to fit the data.  However, this modelling approach would result in $m(p+1)$
unknown coefficients to estimate, which is too many, because $m$ is
usually of the magnitude of hundreds, or even more, in practice.  This
modelling also ignores the similarity which may exist among some individuals.
Such similarity may have very important practical meaning, and could lead to
some important findings in practice.  In addition to that, statistically
speaking, the modelling, like (\ref{eq1}) without any conditions imposed, would
also pay a price on variance side of the estimators resulted because the
available information is not used up.

In order to explore the individual attributes and account for the similarity
among some individuals at the same time, \cite{ke15} proposed a penalised
likelihood/least squares based approach to pursue the homogeneity in the
linear models, i.e. (\ref{eq1}), used for panel data analysis, under the
framework of treating homogeneity as a kind of sparsity.   Regression under
homogeneity condition has also been studied by quite a few recent works, e.g.
 \cite{tibshirani05,friedman07,bondell2008simultaneous,jiangqian13}
, and the references therein.  Like \cite{ke15},
the methods in these works are all based on penalised likelihood/least
squares.  \cite{ke16} took a different approach, they formulated the
homogeneity pursuit problem as a problem of change point detection and applied
the binary segmentation approach to identify the homogeneity in the
linear models with interactive effects.

The existing literature about homogeneity pursuit mainly focuses on the linear
models.  It is well known that the linearity condition may not hold for many
datasets, and the exploration of linear relationship is not sufficient in many
cases.  As a consequence, the semiparametric modelling is becoming more and
more useful in panel data analysis.  Among various semiparametric models, the
single index models
have many advantages, and are a very successful tool in data analysis, see
 \cite{hardle1989investigating,carroll97,yuruppert02,zhuxue06,xia08,penghuang11,zhumiaopeng12,guo2016dynamic},
 and the
reference therein.  In this paper, we are going to investigate the homogeneity
pursuit in the single index models used for panel data analysis.  The detailed
definition of the models we are going to address in this paper is given in
Section \ref{mo0}

\subsection{The single index models with homogeneity structure}
\label{mo0}

Let $y_{it}$ and ${\bf X}_{it}$, a $(p+1)$-dimensional vector, be
respectively the $t$th observations of the response variable and covariate of
the $i$th individual, $i=1, \ \cdots, \ m; \ \  t=1, \ \cdots, \ T$. We
consider the models
\begin{equation}
y_{it}
=
g_i({\bf X}_{it}^{\rm T}\mbox{\boldmath $\beta$}_i) + \epsilon_{it},
\quad
i=1, \ \cdots, \ m; \ \  t=1, \ \cdots, \ T,
\label{model}
\end{equation}
where
\begin{equation}
g_i(\cdot)
=
\left\{
\begin{array}{ll}
\ g_{(1)}(\cdot) & \mbox{ when }  i \in G_{1,1},
\\
\ g_{(2)}(\cdot) & \mbox{ when }  i \in G_{1,2},
\\
\quad \vdots & \quad \vdots
\\
\ g_{(H_1)}(\cdot) & \mbox{ when }  i \in G_{1,H_1},
\end{array}
\right.
\quad
\beta_{ij}
=
\left\{
\begin{array}{ll}
\ \beta_{(1)} & \mbox{ when }  (i, j) \in G_{2,1},
\\
\ \beta_{(2)} & \mbox{ when }  (i, j) \in G_{2,2},
\\
\quad \vdots & \quad \vdots
\\
\ \beta_{(H_2)} & \mbox{ when }  (i, j) \in G_{2,H_2},
\end{array}
\right.
\label{con}
\end{equation}
$\mathbb{G}_1=\{G_{1,k}: \ k = 1, \ \cdots, \ H_1\}$ is a partition of set
$\{1, \ \cdots, \ m\}$, $\mathbb{G}_2=\{G_{2,k}: \ k = 1, \ \cdots, \ H_2\}$ is a partition
of set $\{(i,j): \ i = 1, \ \cdots, \ m; \ j=1, \ \cdots, \ p\}$, $\beta_{ij}$
is the $(j+1)$th component of $\mbox{\boldmath $\beta$}_i$,
and
$$
E(\epsilon_{it}| {\bf X}_{it}) = 0,
\quad
\mbox{var}(\epsilon_{it}| {\bf X}_{it}) = \sigma^2.
$$
The condition (\ref{con}) is the homogeneity structure of the standard single
index models for panel data analysis.  $\{G_{1,k}: \ k = 1, \ \cdots, \ H_1\}$
and $\{G_{2,k}: \ k = 1, \ \cdots, \ H_2\}$ are unknown partitions.  $H_1$ and
$H_2$ are unknown integers, $H_1$ is much smaller than $m$, $H_2$ is much
smaller than $mp$.  $g_{(k)}(\cdot)$, $k=1, \ \cdots, \ H_1$, are unknown
functions to be estimated, and $\beta_{(k)}$, $k=1, \ \cdots, \ H_2$, are
unknown parameters to be estimated.

Let $\beta_{i0}=1$ be the first component of $\mbox{\boldmath $\beta$}_i$.  In the
literature, the most commonly used identification condition for the single
index models is $\|\mbox{\boldmath $\beta$}_i\|=1$ and $\beta_{i0}>0$, or $\beta_{i0}=1$.  We
choose the latter in this paper.

The models (\ref{model}) together with (\ref{con}) show that the homogeneity
pursuit in the single index models for panel data analysis is even more
important than that in the linear models, this is because we would have to
estimate $m$ unknown functions and $mp$ unknown parameters in order to explore
the individual attributes, if the homogeneity pursuit is not conducted.
However, if the homogeneity pursuit is conducted, we only need to estimate
$H_1$, much smaller than $m$, unknown functions and $H_2$, much smaller than
$mp$, unknown parameters when the homogeneity exists.  Even without taking into
account the benefit resulted from the homogeneity pursuit for the parametric
part of the models, just for the part of unknown functions alone, to estimate
much fewer functions would make a big difference in the obtained estimators, in
terms of the stability of the estimators.

The rest of the paper is organized as follows.  We begin in Section \ref{est}
with a description of the proposed estimation procedure which is embedded with
a binary segmentation based homogeneity pursuit.  The asymptotic properties of
the proposed estimators are presented in Section \ref{asy}.  The performance of
the proposed estimation procedure and homogeneity pursuit method, when sample
size is finite, are assessed by simulation studies in Section \ref{simu}.  In
Section \ref{re}, applying the single index models (\ref{model}) together with
the homogeneity structure (\ref{con}) to the 49 Industry Portfolios data set,
which can be freely downloaded from Kenneth French's website
\begin{center}
{\tt
http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data\_library.html},
\end{center}
and the UK climate data, which can be freely downloaded from
\begin{center}
{\tt http://www.metoffice.gov.uk/public/weather/climate-historic},
\end{center}
we will show the advantages of the proposed statistical methodology.  We leave
all technical proofs of the asymptotic properties in the Appendix.

\section{Estimation procedure}
\label{est}

\subsection{Estimation method}
\label{est1}

Our approach to deal with the unknown functions $g_i(\cdot)$,
$i=1, \ \cdots, \ m$, in (\ref{model}) is based on the B-Spline.  To achieve
the best result
for the homogeneity pursuit, we have to decompose all $g_i(\cdot)$s by the same
B-Spline basis, ${\bf B}(\cdot) = (B_1(\cdot), \ \cdots, \ B_K(\cdot))^{{\mbox{\rm\tiny T}}}$.

For each $i$, $i = 1, \ \cdots, \ m$, let $\tilde{\mbox{\boldmath $\beta$}}_i$ be the estimate
of $\mbox{\boldmath $\beta$}_i$ obtained, based on the observations for the $i$th individual,
by a standard estimation procedure for the single index models, e.g. the
method in \cite{yuruppert02} or in \cite{hardle1989investigating}, and
$$
a
=
\min\limits_{1 \leq i \leq m} \min\limits_{1 \leq t \leq T}
{\bf X}_{it}^{{\mbox{\rm\tiny T}}} \tilde{\mbox{\boldmath $\beta$}}_i,
\quad
b
=
\max\limits_{1 \leq i \leq m} \max\limits_{1 \leq t \leq T}
{\bf X}_{it}^{{\mbox{\rm\tiny T}}} \tilde{\mbox{\boldmath $\beta$}}_i.
$$
We use the B-Spline basis of order $s$ in this paper, and the basis,
${\bf B}(\cdot)$, is formed by the equally spaced knots, $\tau_k$,
$k=0, \ \cdots, \ K-s+1$, on the interval $[a, \ b]$, with $\tau_0 = a$ and
$\tau_{K-s+1} = b$.  Based on the basis ${\bf B}(\cdot)$, $g_i(\cdot)$ can be
decomposed as
\begin{equation}
g_i(\cdot) \approx {\bf B}(\cdot)^{{\mbox{\rm\tiny T}}} \mbox{\boldmath $\theta$}_i,
\label{app}
\end{equation}
where $\mbox{\boldmath $\theta$}_i = (\theta_{i1}, \ \cdots, \ \theta_{iK})^{{\mbox{\rm\tiny T}}}$.  So, to get
the estimator of $g_i(\cdot)$, we only need to get the estimator of $\mbox{\boldmath $\theta$}_i$.

Our estimation procedure for $\mbox{\boldmath $\theta$}_i$ and $\mbox{\boldmath $\beta$}_i$, $i=1, \ \cdots, \ m$,
consists of three stages: in the first stage, for each $i$, we estimate
$\mbox{\boldmath $\theta$}_i$ and $\mbox{\boldmath $\beta$}_i$ only based
on the observations for the $i$th individual, and treat the obtained
estimators as initial estimators; we identify, in the second stage, the
homogeneity structure in the $\mbox{\boldmath $\theta$}_i$s and $\mbox{\boldmath $\beta$}_i$s based on the initial
estimators obtained in the first stage; in the final stage, we estimate the
$\mbox{\boldmath $\theta$}_i$s and $\mbox{\boldmath $\beta$}_i$s under the identified homogeneity structure.

We now present the details of the estimation procedure.
\begin{enumerate}
\item[{\bf Stage 1}] ({\it Initial Estimation}).  Let
$\overline\mbox{\boldmath $\beta$}_i=(\beta_{i1}, \ \cdots, \ \beta_{ip})^{\rm T}$, which
is $\mbox{\boldmath $\beta$}_i$ with the first component, which is always $1$, being dropped.
For each $i$, based on the observations for the $i$th individual, approximating
$g_i(\cdot)$ by its decomposition (\ref{app}) and applying the least squares
estimation method, we have the following objective function
\begin{equation}
\sum_{t=1}^T
\Big(
y_{it} - {\bf B}^{\rm T}({\bf X}_{it}^{\rm T}\mbox{\boldmath $\beta$}_i)\mbox{\boldmath $\theta$}_i
\Big)^2.
\label{eqn:min1}
\end{equation}
Minimise (\ref{eqn:min1}) with respect to
$(\overline\mbox{\boldmath $\beta$}_i^{{\mbox{\rm\tiny T}}}, \ \mbox{\boldmath $\theta$}_i^{{\mbox{\rm\tiny T}}})$, and denote the resulting minimiser by
$(\tilde{\mbox{\boldmath $\beta$}}_i^{{\mbox{\rm\tiny T}}}, \ \tilde{\mbox{\boldmath $\theta$}}_i^{{\mbox{\rm\tiny T}}})$.  We will show how to
conduct the minimisation in Section \ref{est2}.

\item[{\bf Stage 2}] ({\it Homogeneity Pursuit}).  Let $\tilde{\beta}_{ij}$ be
the $j$th component of
$\tilde{\mbox{\boldmath $\beta$}}_i$, we sort $\tilde{\beta}_{ij}$, $i=1, \ \cdots, \ m$,
$j=1, \ \cdots, \ p$, in ascending order, and denote them by
$$
b_{(1)}
\leq
\cdots
\leq
b_{(mp)}.
$$
We use $R_{ij}$ to denote the rank of $\tilde{\beta}_{ij}$.  Identifying the
homogeneity among $\tilde{\beta}_{ij}$, $i=1, \ \cdots, \ m$,
$j=1, \ \cdots, \ p$, is equivalent to detecting the change points among
$b_{(l)}$, $l=1, \ \cdots, \ mp$.  To this end, we apply the Binary
Segmentation algorithm as follows.

For any $1 \leq i <  j \leq mp$, let
$$
\Delta_{ij}(\kappa)
=
\sqrt{\frac{(j-\kappa)(\kappa-i+1)}{j-i+1}}
\left|
\frac{\sum_{l=\kappa+1}^j b_{(l)}}{j-\kappa}
-
\frac{\sum_{l=i}^{\kappa} b_{(l)}}{\kappa-i+1}
\right|.
$$
Given a threshold $\delta$, the Binary Segmentation algorithm to detect the
change points works as follows
\begin{enumerate}
\item[(1)] Find $\hat{k}_1$ such that
$$
\Delta_{1,mp}(\hat{k}_1)
=
\max\limits_{1 \leq \kappa < mp} \Delta_{1,mp}(\kappa).
$$
If $\Delta_{1,mp}(\hat{k}_1) \leq \delta$, there is no change point among
$b_{(l)}$, $l=1, \ \cdots, \ mp$, and the process of detection ends.
Otherwise, add $\hat{k}_1$ to the set of change points and divide the region
$\{\kappa:  \ 1 \leq \kappa \leq mp \}$ into two subregions:
$\{\kappa:  \ 1 \leq \kappa \leq \hat{k}_1 \}$ and
$\{\kappa:  \ \hat{k}_1+1 \leq \kappa \leq mp \}$.

\smallskip
\item[(2)] Detect the change points in the two subregions obtained in (1),
respectively.  Let us deal with the region $\{\kappa:  \ 1 \leq \kappa \leq
\hat{k}_1 \}$ first.  Find $\hat{k}_2$ such that
$$
\Delta_{1,\hat{k}_1}(\hat{k}_2)
=
\max\limits_{1 \leq \kappa < \hat{k}_1} \Delta_{1,\hat{k}_1}(\kappa).
$$
If $\Delta_{1,\hat{k}_1}(\hat{k}_2) \leq \delta$, there is no change point in
the region $\{\kappa:  \ 1 \leq \kappa \leq \hat{k}_1 \}$.  Otherwise,
add $\hat{k}_2$ to the set of change points and divide the region
$\{\kappa:  \ 1 \leq \kappa \leq \hat{k}_1 \}$ into two subregions:
$\{\kappa:  \ 1 \leq \kappa \leq \hat{k}_2 \}$ and
$\{\kappa:  \ \hat{k}_2+1 \leq \kappa \leq \hat{k}_1 \}$.  For the region
$\{\kappa:  \ \hat{k}_1+1 \leq \kappa \leq mp \}$, we find $\hat{k}_3$ such
that
$$
\Delta_{\hat{k}_1+1,mp}(\hat{k}_3)
=
\max\limits_{\hat{k}_1+1 \leq \kappa < mp} \Delta_{\hat{k}_1+1,mp}(\kappa).
$$
If $\Delta_{\hat{k}_1+1,mp}(\hat{k}_3) \leq \delta$, there is no change point in the region
$\{\kappa:  \ \hat{k}_1+1 \leq \kappa \leq mp \}$.  Otherwise, add $\hat{k}_3$
to the set of change points and divide the region $\{\kappa:  \ \hat{k}_1+1
\leq \kappa \leq mp \}$ into two subregions:
$\{\kappa:  \ \hat{k}_1+1 \leq \kappa \leq \hat{k}_3 \}$ and
$\{\kappa:  \ \hat{k}_3+1 \leq \kappa \leq mp \}$.
\smallskip
\item[(3)] For each subregion obtained in (2), we do exactly the same as that
for the subregion $\{\kappa:  \ 1 \leq \kappa \leq \hat{k}_1 \}$ or
$\{\kappa:  \ \hat{k}_1+1 \leq \kappa \leq mp \}$ in (2), and
keep doing so until there is no subregion containing any change point.
\end{enumerate}

We sort the estimated change point locations in ascending order and denote them
by
$$
\hat{k}_{(1)} < \hat{k}_{(2)} < \cdots < \hat{k}_{(\hat{H}_{-1})},
$$
where $\hat{H}_{-1}$ is the number of change points detected.  In addition, we
denote $\hat{k}_{(0)} = 0$, $\hat{H}_2 = \hat{H}_{-1} + 1$, and
$\hat{k}_{(\hat{H}_2)} = mp$.

We use $\hat{H}_2$ to estimate $H_2$.  Let
$$
\hat{G}_{2, s}
=
\{(i, j): \ \hat{k}_{(s-1)} < R_{ij} \leq \hat{k}_{(s)}\},
\quad
1 \leq s \leq \hat{H}_2,
$$
we use $\left\{\hat{G}_{2, s}: \ 1 \leq s \leq \hat{H}_2 \right\}$
to estimate the partition $\{G_{2,s}: \ 1 \leq s \leq H_2\}$.  We consider all
the $\beta_{ij}$s with the subscript $(i,j)$ in the same member of the
estimated partition having the same value.

Let $\tilde{\theta}_{ij}$ be the $j$th component of $\tilde{\mbox{\boldmath $\theta$}}_i$.  Doing
exactly the same to $\tilde{\theta}_{ij}$, $i=1, \ \cdots, \ m$, $j=1, \
\cdots, \ K$, we get a partition $\{\hat G_{1,1}, \ \cdots, \ \hat G_{1,\hat{H}_1}\}$ of
$\{(i, j): \ i=1, \ \cdots, \ m; \ j=1, \ \cdots, \ K\}$.  We consider all the
$\theta_{ij}$s with subscript $(i,j)$ in the same member of the estimated
partition having the same value.

\item[{\bf Stage 3}] ({\it Final Estimation}).  Let
$L(\eta_1, \ \cdots, \ \eta_{\hat{H}_2}, \ \xi_1, \ \cdots, \ \xi_{\hat H_1})$
be
\begin{equation}
\label{eqn:min2}
\sum\limits_{i=1}^m \sum_{t=1}^T
\Big(
y_{it} - {\bf B}^{\rm T}({\bf X}_{it}^{\rm T}\mbox{\boldmath $\beta$}_i)\mbox{\boldmath $\theta$}_i
\Big)^2.
\end{equation}
with $\beta_{ij}$, $i=1, \ \cdots, \ m$, $j=1, \ \cdots, \ p$, being replaced
by $\eta_k$ if $(i, j) \in \hat{G}_{2, k}$, and $\theta_{ij}$, $i=1, \ \cdots,
\ m$, $j=1, \ \cdots, \ K$, being replaced by $\xi_s$ if
$(i, j) \in \hat G_{1,s}$.  Let
$(\hat{\eta}_1, \ \cdots, \ \hat{\eta}_{\hat{H}_2}, \ \hat{\xi}_1, \ \cdots, \
\hat{\xi}_{\hat H_1})$ minimise
$L(\eta_1, \ \cdots, \ \eta_{\hat{H}_2}, \ \xi_1, \ \cdots, \ \xi_{\hat H_1})$.
The final estimator $\hat{\beta}_{ij}$ of $\beta_{ij}$ is $\hat{\eta}_k$ if
$(i, j) \in \hat{G}_{2, k}$, and the final estimator $\hat{\theta}_{ij}$ of
$\theta_{ij}$ is $\hat{\xi}_s$ if $(i, j) \in \hat G_{1,s}$.  Once we have the
estimator $\hat{\theta}_{ij}$, the estimator $\hat{g}_i(\cdot)$ of $g_i(\cdot)$
is taken to be ${\bf B}(\cdot)^{{\mbox{\rm\tiny T}}} \hat{\mbox{\boldmath $\theta$}}_i$.

\end{enumerate}

\bigskip

\noindent
\begin{rmk}  When dealing with the unknown functions $g_i(\cdot)$,
$i=1, \ \cdots, \ m$, in the estimation procedure, instead of treating each
unknown function as a single undivided unit to conduct homogeneity pursuit, we
work on the coefficients of its B-Spline decomposition.  This is because there
may still be some kind of homogeneity between two functions even if they are
different.  For example, for two different functions, it could be the case that
some coefficients of the B-Spline decomposition of one function are the same as
some coefficients of the B-Spline decomposition of another one.  If we treat
each unknown function as a single undivided unit to conduct homogeneity
pursuit, we would not identify or use this kind of homogeneity, which would
make our final estimators not as efficient as they should.
\end{rmk}

\subsection{Computational algorithm}
\label{est2}

In the estimation procedure described in Section \ref{est1}, the minimiser of
(\ref{eqn:min1}) does not have a closed form, neither does the minimiser of
$L(\eta_1, \ \cdots, \ \eta_{\hat{H}_2}, \ \xi_1, \ \cdots, \ \xi_{\hat H_1})$.  To
conduct the minimisation of either of the two objective functions, we appeal to
the standard NLS algorithm, and use the {\tt nlsLM of minpack.lm} package in R
to implement it.   One
can also use other NLS software, for example, the NLS routine {\tt lsqnonlin()}
from MATLAB and {\tt PROC NLIN} from SAS.  To use the {\tt nlsLM of minpack.lm}
package in R, we first need to find an initial value.  The initial value for
minimising (\ref{eqn:min1}) can be obtained as follows:
\begin{enumerate}
\item[(1)] Apply the standard least squares estimation for the linear models to
$(y_{it}, \ {\bf X}_{it})$, $t = 1, \ \cdots, \ T$, and denote the resulting
estimator by $\check{\mbox{\boldmath $\beta$}}_i$, the initial value for $\mbox{\boldmath $\beta$}_i$ is taken to
be $\mbox{\boldmath $\beta$}_i^{(0)} = \check{\beta}_{i0}^{-1} \check{\mbox{\boldmath $\beta$}}_i$,
$\check{\beta}_{i0}$ is the first component of $\check{\mbox{\boldmath $\beta$}}_i$.

\item[(2)] Substitute $\mbox{\boldmath $\beta$}_i^{(0)}$ for $\mbox{\boldmath $\beta$}_i$ in (\ref{eqn:min1}),
then minimise (\ref{eqn:min1}) with respect to $\mbox{\boldmath $\theta$}_i$, the minimiser
$\mbox{\boldmath $\theta$}_i^{(0)}$ is the initial value of $\mbox{\boldmath $\theta$}_i$.
\end{enumerate}
Once we have $\mbox{\boldmath $\beta$}_i^{(0)}$ and $\mbox{\boldmath $\theta$}_i^{(0)}$, the minimiser of
(\ref{eqn:min1}) can be obtained by the {\tt nlsLM of minpack.lm} package in R
straightforwardly.

For any set $A$, let $|A|$ be the number of elements in $A$.  The initial
value for minimising $L(\eta_1, \ \cdots, \ \eta_{\hat{H}_2}, \ \xi_1, \
\cdots, \ \xi_{\hat H_1})$ can
be obtained through the initial estimates of $\mbox{\boldmath $\beta$}_i$ and $\mbox{\boldmath $\theta$}_i$,
obtained in Stage 1 of the estimation procedure in Section \ref{est1}, as
follows:
$$
\eta_s^{(0)}
=
\left(|\hat{G}_{2,s}| \right)^{-1}
\sum\limits_{(i,j) \in \hat{G}_{2,s}} \tilde{\beta}_{ij},
\quad
s = 1, \ \cdots, \ \hat{H}_2
$$
and
$$
\xi_s^{(0)}
=
\left(|\hat G_{1,s}| \right)^{-1}
\sum\limits_{(i,j) \in \hat G_{1,s}} \tilde{\theta}_{ij},
\quad
s = 1, \ \cdots, \ \hat H_1.
$$
Once we have the initial value $(\eta_1^{(0)}, \ \cdots, \
\eta_{\hat{H}_2}^{(0)}, \ \xi_1^{(0)}, \ \cdots, \ \xi_{\hat H_1}^{(0)})$, we
can have the minimiser of $L(\eta_1, \ \cdots, \ \eta_{\hat{H}_2}, \
\xi_1, \ \cdots, \ \xi_{\hat H_1})$ by using the {\tt nlsLM of minpack.lm}
package in R straightforwardly.

\subsection{Selection of tuning parameters}
\label{tune}

The threshold $\delta$ in the Stage 2 of the proposed estimation procedure,
described in Section \ref{est1}, plays a key role for the success of the
homogeneity pursuit.  As far as the implementation of the homogeneity pursuit
is concerned, the selection of $\delta$ is equivalent to the selection of
$\hat H_1$ and $\hat H_2$, and to select an integer is easier, therefore,
in this section, instead of selecting $\delta$, we develop a cross-validation
procedure to select the two tuning parameters, $\hat H_1$ and $\hat H_2$.

For the single index model (\ref{model}) where $\mbox{\bf X}_{it}$'s are independent
across $t=1, \ \cdots, \ T,$  we implement a $L$-fold cross validation
approach.  In particular, for a given pair $\{H_1, \ H_2\}$, we remove
$1/L$th of the observed time points for
$\{(y_{it}, \ {\bf X}_{it}), \ i=1, \ \cdots, \ m, \ t=1, \ \cdots, \ T\}$
as a validation set, estimate the single index model (\ref{model}) with
identified homogeneity structure on the remaining data, compute the squared
error between $y_{it}$ and fitted values
$\hat g_i({\bf X}_{it}^{{\mbox{\rm\tiny T}}}\hat\mbox{\boldmath $\beta$}_i)
={\bf B}^{\rm T}({\bf X}_{it}^{\rm T}\hat\mbox{\boldmath $\beta$}_i)\hat\mbox{\boldmath $\theta$}_i$, on the validation set,
and repeat this procedure $L$ times to calculate the cross-validated mean
squared error and its corresponding standard error.  We search over a grid of
$\{H_1, \ H_2\}$ values and apply the one-standard-error rule to choose the
smallest model for which the estimated cross-validated error is within one
standard error of the lowest point on the error surface.  The rationale here is
that if a set of models appear to be more or less equally good, then we might
tend to choose the simplest model.  Across the candidate pairs,
$\{\hat H_1, \ \hat H_2\},$ whose corresponding errors are within this
deviation, one can choose the smallest $\hat H_1$ after selecting the smallest
$\hat H_2$ or switch the selection order or select the smallest value of
$\hat H_1+\hat H_2$, we take the first approach since it produces better model
selection consistency in our numerical experiments.  A similar
one-standard-deviation-rule technique has been adopted to choose the
regularisation parameter with a smaller model size for the lasso problems
\cite[]{james13}.

When $\mbox{\bf X}_{it}$'s are time dependent panel data, we implement a rolling
procedure to perform cross-validation for time series.  More specifically,
for each $r=L, \ L-1, \ \cdots, \ 1,$ we rollingly treat
$\{(y_{it}, \ {\bf X}_{it}), \ i=1, \ \cdots, \ m, \ t=1, \ \cdots, \ T-r\}$ as
training observations and $\{(y_{i,T-r+1}, \ {\bf X}_{i,T-r+1}), \ i=1, \ \cdots, \
m\}$ as validation set, calculate the squared error between each $y_{it}$ and
its fitted value.  Finally, we apply the one-standard-deviation-rule on the
lowest cross-validated mean squared error and choose $\hat H_2$ and $\hat H_1.$

In the cross-validation procedure when we need make predictions for validation
set, the domain in ${\bf B}({\bf X}_{it}^{{\mbox{\rm\tiny T}}}\hat \mbox{\boldmath $\beta$}_i)$ for traning data set might
not cover that for validation set.  We adopt the idea in \cite{wangyang09} by
mapping ${\bf X}^{{\mbox{\rm\tiny T}}}_{it}\mbox{\boldmath $\beta$}_i$ to
$F_i({\bf X}^{{\mbox{\rm\tiny T}}}_{it}\mbox{\boldmath $\beta$}_i) \sim \text{Unif}[0, \ 1]$,
where $F_i$ is the distribution function of ${\bf X}^{{\mbox{\rm\tiny T}}}_{it}\mbox{\boldmath $\beta$}_i$.  We then
implement the estimation procedure described in Section~\ref{est1} by
decomposing
$g_i({\bf X}^{{\mbox{\rm\tiny T}}}_{it}\mbox{\boldmath $\beta$}_i)
\approx
{\bf B}^{{\mbox{\rm\tiny T}}}\big(F_i({\bf X}^{{\mbox{\rm\tiny T}}}_{it}\mbox{\boldmath $\beta$}_i)\big) \mbox{\boldmath $\theta$}_i$.
The proposed approach is thus able to make predictions and, as demonstrated by
some numerical studies, provides very similar sample performance in terms of
estimation accuracy.

\subsection{Post-processing step}
\label{post}
We equip the Binary-Segmentation-algorithm-based homogeneity pursuit with an
additional step aimed to enhance the accuracy of detected change-points
locations through a fine-scale search. To be specific, at each change-point,
we re-calculate  $\Delta_{ij}(\kappa)$ over the interval between two adjacent
change-points and identify the new change-point location to replace the old
one. We perform this post-processing procedure by iteratively cycling
through all neighbouring change-points and fine-tuning the change-points
locations. This procedure is terminated when the set of change-points does not
change. Our numerical experiments show that this extra post-processing step
apparently improve the accuracy of each estimated change-point location and
hence the identified homogeneity structure for model~\eqref{model}.

\section{Asymptotic properties}
\label{asy}

In this section, we are going to investigate the asymptotic behaviour of the
estimators obtained by the proposed estimation procedure, which we call
correct-fitting, and compare with the estimators obtained without homogeneity
pursuit, which is the initial estimators obtained in the Stage 1 in the
proposed estimation procedure, we call it over-fitting, and the estimators
obtained under the assumption that all individuals share the same index
(namely, $\mbox{\boldmath $\beta$}_1 = \cdots = \mbox{\boldmath $\beta$}_m$), which we call under-fitting.  The
asymptotic theory presented in this section is in the sense that
$T\longrightarrow \infty$, and $m$, $p$ are all possibly diverging to infinity
but $H_1$, $H_2$ are fixed.  This agrees with many applications in which $H_1$
and $H_2$ are expected to be small and thus significant reduction of unknown
parameters can be achieved by clustering the parameters.  To make the
presentation neat, we state the asymptotic theorems in this section and leave
all technical proofs in the Appendix.

Let $\epsilon_{t}=(\epsilon_{1t}, \ \cdots, \ \epsilon_{mt})^{\rm T}$,
$$
\overline{\bf X}_{it}=(X_{it,1}, \ \cdots, \ X_{it,p})^{\rm T},
\quad
y_t=(y_{1t}, \ \cdots, \ y_{mt})^{\rm T},
\quad
{\bf X}_{t}=({\bf X}_{1t}^{\rm T}, \ \cdots, \ {\bf X}_{mt}^{\rm T})^{\rm T}.
$$
In this paper, we assume $(y_{t}, \ {\bf X}_{t}, \ \epsilon_t)$ are stationary
with $\alpha(l)\le \rho^l$ for some $\rho<1$, and $\epsilon_t$ is
independent of ${\bf X}_t$.  Note that unlike \cite{vogt15}, we do not need to
assume independence or stationarity of variables cross $i$.

We start with the asymptotic properties of the estimators obtained without
homogeneity pursuit. The convergence rate of the estimator $\tilde{\beta}_{ij}$ is of order $T^{-1/2}$, and the convergence rate of the estimator
$\tilde{g}_i(u)$ is of order $T^{-2/5}$, which is as expected as we assumed the
functions are twice differentiable.

\bigskip

\begin{thm}\label{thm:new1}
({\it Over-fitting case}).  For any $i$, $i=1, \ \cdots, \ m$, and
$1\le j\le p$, under the conditions (C1)-(C4) and (C5') in the Appendix, we
have
$$
T^{1/2} (\widetilde\varepsilon_{ij}^{\rm T}\widetilde\mbox{\boldmath $\Theta$}_2\widetilde\varepsilon_{ij})^{-1/2}
\left(\tilde{\beta}_{ij} - \beta_{ij} \right)
\stackrel{D}{\longrightarrow} N(0,\ 1)
$$
and
$$
T^{2/5}(\widetilde{\bf b}_i^{\rm T}(u)\widetilde\mbox{\boldmath $\Theta$}_1\widetilde{\bf b}_i(u))^{-1/2}
\left(\widetilde{g}_i(u) - g_{i}(u)-r_i(u) \right)
\stackrel{D}{\longrightarrow} N(0,\ 1),
$$
where $\widetilde\varepsilon_{ij}$ and $\widetilde{\bf b}_i(u)$ are unit vectors,  $\widetilde\mbox{\boldmath $\Theta$}_1$,
$\widetilde\mbox{\boldmath $\Theta$}_2$ are matrices with eigenvalues bounded and bounded away from
zero, all these quantities are defined in the proof in the
Appendix \ref{sec:conv4}.  The bias term
$r_i(u)=g_{i}(u)-{\bf B}^{\rm T}(u)\mbox{\boldmath $\theta$}_{0i}$ satisfies $|r_i(u)|\le CK^{-2}$,
where $\mbox{\boldmath $\theta$}_{0i}$ is the vector of spline coefficients used to approximate
$g_i$ as defined in Appendix \ref{sec:ass}
parameter of $g_i$ as defined in assumption  (C3).
\end{thm}

\bigskip

Let $m_i$ be the size of $G_{1,h}$ that contains $i$, and $m_{ij}$ be the size
of $G_{2,h}$ that contains $\beta_{ij}$.  To make the statement about the
correct-fitting case cleaner, we assume that all $m_i$ are of the same order
and all $m_{ij}$ are of the same order ($\max_{i,j} m_{ij}/\min_{i,j} m_{ij}$
and $\max_i m_i/\min_i m_i$ are bounded) in the following theorem, which shows
in particular that the convergence rate of the estimator $\hat{\beta}_{ij}$ is
of order $(mpT)^{-1/2}$, and the convergence rate of the estimator
$\hat{g}_i(u)$ is of order $(mT)^{-2/5}$.

\bigskip

\begin{thm}\label{thm:new2} ({\it Correct-fitting case}).  For any $i$,
$i=1, \ \cdots, \ m$, and $1\le j\le p$, under the conditions (C1)-(C6) in the
Appendix, we have
$$
(mpT)^{1/2} (\varepsilon_{ij}^{\rm T}\mbox{\boldmath $\Theta$}_2\varepsilon_{ij})^{-1/2}
\left(\hat{\beta}_{ij} - \beta_{ij} \right)
\stackrel{D}{\longrightarrow} N(0, \ 1)
$$
and
$$
(mT)^{2/5}({\bf b}_i^{\rm T}(u)\mbox{\boldmath $\Theta$}_1{\bf b}_i(u))^{-1/2}
\left(\hat{g}_i(u) - g_i(u)-r_i(u) \right)
\stackrel{D}{\longrightarrow} N(0, \ 1),
$$
where $\varepsilon_{ij}$ and ${\bf b}_i(u)$ are unit vectors, $\mbox{\boldmath $\Theta$}_1$, $\mbox{\boldmath $\Theta$}_2$ are
matrices with eigenvalues bounded and bounded away from zero, all these
quantities are defined in the proof in the Appendix \ref{sec:conv4}.
\end{thm}

\bigskip

Finally, for the under-fitting case, let $\check{\mbox{\boldmath $\beta$}}_i$ and
$\check{g}_i(\cdot)$ be the estimators of $\mbox{\boldmath $\beta$}_i$ and $g_i(\cdot)$ obtained
under the assumption that all individuals share the same unknown parameters.

\bigskip

\begin{thm}\label{thm:new3} ({\it Under-fitting case}).  Suppose the
$\mbox{\boldmath $\beta$}_i$s are sufficiently separated in the sense that for
$\bar\mbox{\boldmath $\beta$}:=\sum_{i=1}^m\mbox{\boldmath $\beta$}_i/m$,
$$\frac{1}{mp}\sum_{i=1}^m\|\mbox{\boldmath $\beta$}_i-\bar\mbox{\boldmath $\beta$}\|^2\ge c$$
for some $c>0$, then
$$\frac{1}{mp}\sum_{i=1}^m\|\check \mbox{\boldmath $\beta$}_i-\mbox{\boldmath $\beta$}_i\|^2\ge c.$$
Similarly, if
$$\frac{1}{m}\sum_{i=1}^m\int\left |{g}_i(u) - \bar g(u) \right|^2 du \ge c,
$$
where $\bar g(u)=m^{-1}\sum_{i=1}^m g_i(u)$, then
$$\frac{1}{m}\sum_{i=1}^m\int\left |\check{g}_i(u) -  g_i(u) \right|^2 du \ge c.
$$
\end{thm}

\section{Simulation studies}
\label{simu}

In this section, we are going to use a simulated example to demonstrate how
accurate the proposed estimation is.  We will also show much loss it would
inflict if the homogeneity structure is ignored or mistakenly specified as that
all individuals share the same index coefficients or the same link function.

\noindent
{\bf Example}.  We generate a sample from model (\ref{model}) with $p=2$ and
WLOG an even $m,$
$$
g_i(u)
=
\begin{cases}
\sin(\pi u/4) & \text{when } i=1, \ 2, \ \cdots, \ m/2,
\\
\cos(\pi u/4) & \text{when } i=m/2+1, \ \cdots, \ m,
\end{cases}
$$
and
$$
\mbox{\boldmath $\beta$}_i
=
\begin{cases}
(1, \ -1.5\sqrt{0.2}, \ -0.5\sqrt{0.2})^{{\mbox{\rm\tiny T}}} &
\text{when } i=1, \ 3, \ \cdots, \ m-1,
\\
(1, \ 0.5\sqrt{0.2}, \ 1.5\sqrt{0.2})^{{\mbox{\rm\tiny T}}} &
\text{when } i=2, \ 4, \ \cdots, \ m,
\end{cases}
$$
where $\|\mbox{\boldmath $\beta$}_i\|^2=1.5$ for $i=1, \ \cdots, \ m$.  Let ${\bf X}_{it}$ and
$\epsilon_{it}$, $i=1, \ \cdots, \ m$, $t=1, \ \cdots, \ T$ be independently
generated from $\frac{1}{\sqrt{1.5}}N(\mbox{\bf 0}_3, \ {\bf I}_3)$ truncated by
$[-1.343, \ 1.343]^3$ (the range of 5th to 95th quantiles for $N(0, \ 2/3)$)
and $N(0, \ \sigma^2)$, respectively.  Once ${\bf X}_{it}$ and $\epsilon_{it}$ are
generated, $y_{it}$ can be generated through (\ref{model}).

We conduct the simulated example for various $m$s and $T$s with $\sigma=0.2$,
and compare our proposed approach to its potential competitors based on the
following performance metrics:
\begin{enumerate}
\item[(1)]{\it Estimation accuracy}.  For an estimator $\hat{\mbox{\boldmath $\beta$}}_i$ of
$\mbox{\boldmath $\beta$}_i$, we use the mean squared error (MSE), namely
$\mbox{MSE}(\hat{\mbox{\boldmath $\beta$}}_i) = E\left(\|\hat{\mbox{\boldmath $\beta$}}_i - \mbox{\boldmath $\beta$}_i \|^2\right)$, to
assess the estimation
error of $\hat{\mbox{\boldmath $\beta$}}_i$.  Analogously, for an estimator $\hat g_i(\cdot)$ of
$g_i(\cdot)$, its estimation accuracy can be evaluated based on the mean
integrated squared error,
$$
\text{MISE}(\hat g_i)
=
E \left\{\int\left(\hat{g}_i(u) - g_i(u)\right)^2 du \right\}.
$$
To avoid the situation where the performance is dominated by the poor boundary
behaviour, we let the integral domain to be non-boundary region, which is
between the 1st and 99th quantiles of
$\{{\bf X}_{it}^{{\mbox{\rm\tiny T}}}\mbox{\boldmath $\beta$}_i, \ t=1, \ \cdots, \ T\}$.

\item[(2)]{\it Homogeneity structure identification consistency}.  To evaluate
the distance between the detected homogeneity structure and the true one, we
use the normalized mutual information (NMI) \cite[]{ke15}, which measures the
similarity between two partitions.  Suppose
$\mathbb{C}=\{C_1, \ C_2, \ \cdots \}$ and
$\mathbb{D}=\{D_1, \ D_2, \ \cdots \}$ are two partitions of
$\{1, \ \cdots, \ n\}$, the NMI is defined as
$$
\text{NMI}(\mathbb{C}, \ \mathbb{D})
=
\frac{I(\mathbb{C}, \ \mathbb{D})}
{[H(\mathbb{C})+ H(\mathbb{D})]/2},
$$
where
$$
I(\mathbb{C}, \ \mathbb{D})
=
\sum\limits_{k,j}
\big(|C_k \cap D_j|/n\big) \hbox{log}\big(n| C_k \cap D_j|/|C_k||D_j|\big)
$$
and
$$
H(\mathbb{C})
=
-\sum\limits_{k}
\big(|C_k|/n\big) \hbox{log}\big(|C_k|/n\big).
$$
The NMI takes values in $[0, \ 1]$ with larger values indicating higher level
of similarity between two partitions.  For an estimated partition
$\hat{\mathbb{G}}_2=\{\hat G_{2,1}, \ \cdots, \ \hat G_{2, \hat H_2}\}$ of
$\{(i, \ j): \ 1, \ \cdots, \ m, \ j=1, \ \cdots, \ p\}$, obtained in the
Stage 2 of the proposed estimation procedure in Section \ref{est1}, we
calculate $\text{NMI}(\hat{\mathbb{G}}_2, \ \mathbb G_2)$ to assess how close
to the true homogeneity structure in $\beta_{ij}$s the estimated one is.
Similarly, for an estimated partition $\hat{\mathbb{G}}_1$ of
$\{i: \ 1, \ \cdots, \ m\}$, we use
$\text{NMI}(\hat{\mathbb{G}}_1, \ \mathbb G_1)$ to evaluate how close
the estimated homogeneity structure in $g_{i}(\cdot)$s is to the true one.
\end{enumerate}

For each case, we apply either the single index model (\ref{model}) with the
standard estimation procedure, the initial estimation of the proposed
estimation procedure in Section~\ref{est1}, which we call over-fitting (Over),
the single index model (\ref{model}) with the homogeneity structure
(\ref{con}) together with the proposed estimation procedure, which we call
correct-fitting, the single index model (\ref{model}) with all individuals
share the same index vector (namely, $\mbox{\boldmath $\beta$}_1 = \cdots = \mbox{\boldmath $\beta$}_m$),
which we call Under-I, the single index model (\ref{model}) with all
individuals share the same link function (namely,
$g_1(\cdot) = \cdots = g_m(\cdot)$, i.e. $\mbox{\boldmath $\theta$}_1 = \cdots = \mbox{\boldmath $\theta$}_m$),
which we call Under-F, or the single index model (\ref{model}) with all
individuals share both the same index and link function, which we call
Under-I-F, to the simulated data set.

We develop three methods under the correct-fitting case.  The first approach,
named Correct-C, optimises (\ref{eqn:min2}) based on the estimated
{\it componentwise} homogeneity structure in $\beta_{ij}$s and $\theta_{ij}$s,
obtained in the Stage 2 of the proposed estimation procedure in Section
\ref{est1} with the tuning parameters selected through the cross-validation
approach described in Section~\ref{tune}.  The second approach, Correct-V, is
the same as the first approach but optimises (\ref{eqn:min2}) based on the
estimated componentwise homogeneity structure in $\beta_{ij}$s and
{\it vectorwise} homogeneity structure in $\mbox{\boldmath $\theta$}_i$s which can be obtained
through the estimated componentwise homogeneity structure in $\theta_{ij}$s.
The third approach, which we call Correct-NMI, is the same as the second
approach but with the tuning parameters selected to be the one maximising
$\text{NMI}(\hat{\mathbb{G}}_2, \ \mathbb G_2)$ and
$\text{NMI}(\hat{\mathbb{G}}_1, \ \mathbb G_1)$.  In practice without knowing
the true homogeneity structure, one cannot implement Correct-NMI.

Under-I, Under-F and Under-I-F are three kinds of under-fitting.  For Under-I
or Under-F, the homogeneity structure in $\mbox{\boldmath $\theta$}_i$s or $\beta_{ij}$s is
estimated in the same way as that in the proposed estimation
procedure in Section \ref{est1} with the tuning parameters still selected by
the one-standard-deviation-rule cross-validation approach.

We compare over-fitting, correct-fittings and under-fittings to the oracle case
where the true homogeneity structure is used.  The computational algorithms
for the under-fitting and oracle estimators are the same as that for the
correct-fitting, but use either identified or pre-specified homogeneity
structure.  We compare the sample performance of all eight approaches in our
conducted simulation study.

We report the results for estimation errors and NMIs for $\mbox{\boldmath $\beta$}_i$s and
$g_i(\cdot)$s averaged over 100 replicates in Tables~\ref{tab_beta} and
\ref{tab_fun}, respectively.  In terms of estimation error, the overall
estimation accuracy is improved as $m$ and $T$ increase and three
correct-fitting approaches perform very well as reflected in their lower values
of MSEs and MISEs.  Among the three methods, Correct-NMI provides the best
performance even producing very comparable MSEs and MISEs with the oracle
estimator and Correct-C is outperformed by Correct-V in most settings.  This is
somewhat expected, since, unlike Correct-C, which optimises (\ref{eqn:min2})
based on the detected homogeneity structure in $\beta_{ij}$s and
$\theta_{ij}$s, Correct-V separates the final estimation step from the
cross-validation procedure, which is used to identify the homogeneity structure
in $\beta_{ij}$s and $\mbox{\boldmath $\theta$}_{i}$s.   Analogously, Correct-NMI solves a
separate optimisation after detecting the homogeneity structure based on the
largest NMIs.  It is also worth noting that the over-fitting and under-fitting
methods, which either ignores or mistakenly specify the homogeneity structure,
provide much worse results, highlighting the importance of incorporating the
appropriate homogeneity structure.  In terms of selecting the structure of
homogeneity, we observe that three correct-fitting methods produce perfect
identifications of homogeneity structure in $\beta_{ij}$s and Correct-NMI
provides the largest NMI values indicating that it can effectively recover the
true homogeneity structure in $g_i(\cdot)$s.  The performance of Correct-C and
Correct-V deteriorates when $m$ increases, this is intuitively due to the
increased $m$ values and the cross-validation procedure, which tends to choose
a larger number of change points as $m$ increases, resulting in smaller NMI
values for $\hat{\mathbb{G}}_1$.

\begin{table}[htbp]
	\begin{center}
		\caption{{\bf The Average of $\mbox{MSE}(\hat{\mbox{\boldmath $\beta$}}_i)$,
$i=1, \ \cdots, \ m$, and Average NMIs for $\hat{\mathbb{G}}_2$}
\centerline{{\it All entries for MSEs are $10^4$ times their actual values}}}
		\label{tab_beta}
		\vspace{0.3cm}
		\begin{tabular}{cc|ccc|ccc}
			\hline
			\hline
			& $T$ & \multicolumn{3}{c|}{400} & \multicolumn{3}{c}{800}\tabularnewline
			& $m$ & 30 & 60 & 90 & 30 & 60 & 90\tabularnewline
			\hline
			\multirow{8}{*}{MSE}
			& Oracle & 0.380 & 0.205 & 0.133 & 0.214 & 0.128 & 0.101\tabularnewline
			& Correct-C & 0.473 & 0.221 & 0.147 & 0.157 & 0.093 & 0.085 \tabularnewline
			& Correct-V & 0.381 & 0.208 & 0.134 & 0.214 & 0.129 & 0.100 \tabularnewline
			& Correct-NMI & 0.381 & 0.206 & 0.133 & 0.214 & 0.128 & 0.101 \tabularnewline
			& Over & 5.438 & 5.285 & 5.246 & 2.636 & 2.689 & 2.663 \tabularnewline
			& Under-I-F & 4005.6 & 4002.6 & 4002.1 & 4003.4 & 4001.5 & 4001.1\tabularnewline
			& Under-I & 4002.4 & 4001.7 & 4001.3 & 4001.8 & 4001.0 & 4000.6\tabularnewline
			& Under-F & 21.233 & 3.393 & 2.259 & 3.111 & 1.454 & 1.123\tabularnewline
			\hline
			\multirow{8}{*}{NMI}
			& Oracle & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000\tabularnewline
			& Correct-C & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000\tabularnewline
			& Correct-V & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000\tabularnewline
			& Correct-NMI & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000\tabularnewline
			& Over & 0.339 & 0.290 & 0.267 & 0.339 & 0.290 & 0.267\tabularnewline
			& Under-I-F & 0 & 0 & 0 & 0 & 0 & 0\tabularnewline
			& Under-I & 0 & 0 & 0 & 0 & 0 & 0\tabularnewline
			& Under-F & 0.817  & 0.810 & 0.807 & 0.812 & 0.809 & 0.806\tabularnewline
			\hline
			\hline
		\end{tabular}
	\end{center}
\end{table}

\begin{table}[htbp]
	\begin{center}
		\caption{{\bf The Average of $\mbox{MISE}(\hat{g}_i)$, $i=1, \
\cdots, \ m$, and Average NMIs for $\hat{\mathbb{G}}_1$}
\centerline{{\it
where all entries for MISEs are $10^2$ times their actual values
}}}
		\label{tab_fun}
				\vspace{0.3cm}
		\begin{tabular}{cc|ccc|ccc}
			\hline
			\hline
			& $T$ & \multicolumn{3}{c|}{400} & \multicolumn{3}{c}{800}\tabularnewline
			& $m$ & 30 & 60 & 90 & 30 & 60 & 90\tabularnewline
			\hline
			\multirow{8}{*}{MISE}
			& Oracle & 0.260 & 0.249 & 0.243 & 0.251 & 0.241 & 0.236\tabularnewline
			& Correct-C & 0.667 & 0.553 & 0.485 & 0.366 & 0.318 & 0.307 \tabularnewline
			& Correct-V & 0.300 & 0.289 & 0.311 & 0.257 & 0.246 & 0.246\tabularnewline
			& Correct-NMI & 0.266 & 0.254 & 0.249 & 0.252 & 0.241 & 0.238\tabularnewline
			& Over & 0.548 & 0.544 & 0.540 & 0.410 & 0.403 & 0.396\tabularnewline
			& Under-I-F & 87.988 & 87.963 & 87.959 & 88.746 & 88.726 & 88.647\tabularnewline
			& Under-I & 10.099 & 9.853 & 9.952 & 9.770 & 9.719 & 9.711\tabularnewline
			& Under-F & 85.196 & 85.240 & 85.235 & 86.017 & 86.035 & 85.952\tabularnewline
			\hline
			\multirow{8}{*}{NMI}
			& Oracle & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \tabularnewline
			& Correct-C & 0.816 & 0.790 & 0.637 & 0.958 & 0.948 &0.881 \tabularnewline
			& Correct-V & 0.816 & 0.790 & 0.637 & 0.958 & 0.948 & 0.881\tabularnewline
			& Correct-NMI & 0.965 & 0.966 & 0.948 & 0.994 & 0.997 & 0.997\tabularnewline
			& Over & 0.339 & 0.290 & 0.267 & 0.339 & 0.290 & 0.267\tabularnewline
			& Under-I-F & 0 & 0 & 0 & 0 & 0 & 0\tabularnewline
			& Under-I & 0.811 & 0.797 & 0.783 & 0.901 & 0.894 & 0.900\tabularnewline
			& Under-F & 0 & 0 & 0 & 0 & 0 & 0\tabularnewline
			\hline
			\hline
		\end{tabular}
	\end{center}
\end{table}

\section{Real data analysis}
\label{re}

We will illustrate the proposed method with two real data examples in this
section.

\subsection{Industrial Portfolio's return}
\label{port}

We first study the data set about $m=49$ Industrial Portfolios' daily simple
return from 1/8/2015 to 31/12/2015.  This data set can be freely downloaded
from Kenneth French's website

\begin{center}
{\tt http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data\_library.html}
\end{center}

This data set has been analysed in quite a few literature.  For example,
\cite{guo2016dynamic} used this data set to demonstrate the performance of a
newly developed dynamic portfolio allocation.  In this paper, we are going to
explore the homogeneity structure in this data set by our proposed method.

Let $y_{it}$ be the daily simple return of the $i$th portfolio at the $t$th
day, $i=1, \ \cdots, \ m, \ t=1, \ \cdots, \ T$, and
${\bf X}_{it}=(X_{t1}, \ X_{t2}, \ X_{t3})^{{\mbox{\rm\tiny T}}}$ be the observation of the
Fama-French three factors, where $X_{t1}, \ X_{t2}, \ X_{t3}$ respectively
represent the market (Rm-Rf), size (SMB) and value (HML) factors at the $t$th
day.

We apply the single index model (\ref{model}) with the unknown homogeneity
structure (\ref{con}) to fit the data set.  From interpretation point of view,
the homogeneity structure in the unknown link functions, $g_i(\cdot)$s, where
each $g_i(\cdot)$ is treated as a single undivided unit would make much more
sense than the homogeneity structure in the coefficients of the B-Spline
decompositions of $g_i(\cdot)$s.  Therefore, we use the Correct-V, described
in Section~\ref{simu}, to identify the homogeneity structure in $\beta_{ij}$s
or $g_i(\cdot)$s, and estimate the unknown parameters and unknown functions.

In the implementation of the Correct-V, we implement the method in
Section \ref{est1} with the tuning parameters selected by the cross-validation
for time series as described in Section~\ref{tune}.  Specifically, we define
the cross-validated mean squared error
\begin{equation}
\label{cv.err}
\mbox{CV}
=
\frac{1}{m L}\sum\limits_{i=1}^{m} \sum\limits_{t=T-L+1}^{T}
\left(y_{it}-\hat{y}_{it}\right)^2
\end{equation}
where $L=30$.  Note that we here do not apply the one-standard-rule when
performing the cross
validation to select the tuning parameters for identifying the homogeneity
structure, since we have already selected a small enough model with 11 and 2
detected groups in index coefficients and link functions, respectively.
Table~\ref{tab4} provides the identified clustering results for $\beta_{i2}$,
$\beta_{i3}$, $g_{i}(\cdot)$, $i=1, \ \cdots, \ 49$, and
Figure~\ref{funcplot_port} plots the estimated link functions.  We observe a
few apparent patterns.  Firstly, the estimated link functions are very linear
indicating the linear relationship between portfolio returns and Fama-French
three factors, which has been verified by broad empirical studies.  Secondly,
many portfolios belonging to similar industrials were grouped into the same
cluster for the estimated index coefficients, e.g. Hardw, Softw and Agric,
Food, Soda were clustered into Groups~8 and 4 in terms of the estimated
coefficients for factors SMB and HML respectively.

\begin{table}[htbp]
\footnotesize
	\caption{\label{tab4}{\bf
Grouping Results for The Index Coefficients for SMB, HML and
\centerline{Link Functions of 49 Industrial Portfolios}}}
\begin{center}
\begin{tabular}{c|cccccccccc}
	\hline
	\hline
	& Agric & Food & Soda & Beer & Smoke & Toys & Fun & Books & Hshld & Clths\tabularnewline
	SMB & 6 & 5 & 6 & 4 & 2 & 8 & 9 & 7 & 6 & 6\tabularnewline
	HML & 4 & 4 & 4 & 3 & 4 & 5 & 3 & 4 & 5 & 6\tabularnewline
	Function & ii & i & ii & ii & i & i & i & i & i & ii\tabularnewline
	\hline
	& Hth & MedEq & Drugs & Chems & Rubbr & Txtls & BldMt & Cnstr & Steel & FabPr\tabularnewline
	SMB & 9 & 8 & 9 & 7 & 7 & 8 & 8 & 8 & 9 & 10\tabularnewline
	HML & 2 & 2 & 1 & 6 & 4 & 4 & 6 & 6 & 9 & 7\tabularnewline
	Function & i & i & i & i & i & i & i & i & i & ii\tabularnewline
	\hline
	& Mach & ElcEq & Autos & Aero & Ships & Guns & Gold & Mines & Coal & Oil\tabularnewline
	SMB & 8 & 9 & 7 & 6 & 8 & 6 & 10 & 9 & 11 & 9 \tabularnewline
	HML & 8 & 6 & 6 & 6 & 7 & 4 & 10 & 7 & 11 & 9 \tabularnewline
	Function & i & i & i & i & i & i & ii & i & ii & i\tabularnewline
	\hline
	& Util & Telcm & PerSv & BusSv & Hardw & Softw & Chips & LabEq & Paper & Boxes\tabularnewline
	SMB & 4 & 8 & 8 & 7 & 8 & 8 & 7 & 7 & 6 & 6 \tabularnewline
	HML & 5 & 4 & 5 & 4 & 5 & 3 & 5 & 5 & 5 & 6\tabularnewline
	Function & i & i & i & i & i & i & i & i & i & i \tabularnewline
	\hline
	& Trans & Whlsl & Rtail & Meals & Banks & Insur & RIEst & Fin & Other & \tabularnewline
	SMB & 7 & 7 & 7 & 7 & 7 & 6 & 7 & 6 & 7 & \tabularnewline
	HML & 6 & 6 & 4 & 3 & 6 & 5 & 5 & 5 & 6 & \tabularnewline
	Function & i & i & i & i & i & i & i & i & ii & \tabularnewline
	\hline
	\hline
\end{tabular}
\end{center}
\end{table}

\begin{figure*}
	\centering
	\includegraphics[width=16cm,height=16cm]{figure/funcplot1_port.pdf}
	\caption{\label{funcplot_port}{Plots of estimated link functions with
respect to $\{{\bf X}_{it}^{{\mbox{\rm\tiny T}}}\hat \mbox{\boldmath $\beta$}_i, t=1,\dots, T\}$.}}
\end{figure*}

\subsection{UK climate data}
\label{temp}

Our second data set, which is available from the UK Met Office website
\begin{center}
{\tt http://www.metoffice.gov.uk/public/weather/climate-historic},
\end{center}
contains monthly data of the mean daily maximum temperature (TMAX), mean daily
minimum temperature (TMIN), days of air frost (AF), total rainfall (RAIN) and
total sunshine duration (SUN) collected from 37 stations across the UK.  We
first remove the missing values and thus select data during the period of
January 1993 to December 2009 from 16 locations.  We then eliminate the
seasonality and trend effects and standardise the data.  Let $y_{it}$ be the
monthly mean temperature, which can be calculated as (TMAX+TMIN)/2, and
${\bf X}_{it}=(X_{it1}, \ X_{it2}, \ X_{it3})^{{\mbox{\rm\tiny T}}}$ be the observations for AF,
RAIN and SUN, from the $i$th station at the $t$th month, $i=1, \ \cdots, \ 16,
\ t=1, \ \cdots, \ T$,

Like the analysis of the Industrial Portfolio's return data set, we apply the
single index model (\ref{model}) with unknown homogeneity structure (\ref{con})
together with the proposed estimation procedure, Correct-V, to the data set.
Table~\ref{tab6} provides the clustering results for the index coefficients and
link functions, where 4 and 2 groups were selected respectively.
Figure~\ref{funcplot_temp} plots the estimated link functions at 16 stations.

It is very interesting to see, from Table~\ref{tab6}, that Oxford, Hurn,
Eastbourne and Bradford share exactly the same model, which implies the impact
of rainfall or total sunshine duration on monthly mean temperature has exactly
the same pattern in these four areas.  The same finding also appears in the
three areas of Waddington, Sheffield and Heathrow, the two areas of
Ross-On-Wye and Eskdalemuir, and the two areas of Paisley and Leuchars.  If we
only focus on the impact of rainfall on monthly mean temperature, the seven
areas of Waddington, Sheffield, Shawbury, Paisley, Leuchars, Lerwick and
Heathrow would have exactly the same pattern.  Similar finding also appears
for the impact of total sunshine duration on monthly mean temperature.

\begin{table}[htbp]
\footnotesize
	\begin{center}
	\caption{{\bf Grouping Results for The Index Coefficients for RAIN,
SUN and
\centerline{Link Functions at 16 Locations}}}
	\label{tab6}
			\vspace{0.3cm}
		\begin{tabular}{c|cccccccc}
			\hline
			\hline
			& Waddington & Sheffield & Shawbury & Ross-On-Wye & Paisley & Oxford & Leuchars & Lerwick\tabularnewline
			RAIN & 4 & 4 & 4 & 3 & 4 & 3 & 4 & 4 \tabularnewline
			SUN & 1 & 1 & 2 & 2 & 3 & 2 & 3 & 4 \tabularnewline
			Function & i & i & i & ii & i & i & i & i \tabularnewline
			\hline
			& Hurn & Heathrow & Eskdalemuir & Eastbourne & Cambridge & Camborne & Bradford & Armagh \tabularnewline
			RAIN & 3 & 4 & 3 & 3 & 3 & 3 & 3 & 4 \tabularnewline
			SUN & 2 & 1 & 2 & 2 & 1 & 3 & 2 & 2 \tabularnewline
			Function & i & i & ii & i & i & i & i & ii \tabularnewline
			\hline
			\hline
		\end{tabular}
	\end{center}
\end{table}

\begin{figure*}
	\label{fig3}
	\centering
	\includegraphics[width=16cm,height=16cm]{figure/funcplot_temp.pdf}
	\caption{\label{funcplot_temp}{Plots of estimated link functions with
respect to $\{{\bf X}_i^{{\mbox{\rm\tiny T}}}\hat \mbox{\boldmath $\beta$}_{it}, t=1, \dots, T\}$.}}
\end{figure*}

redpajama_set_name: RedPajamaArXiv