% interactnlmsample.tex
% v1.05 - August 2017

\documentclass[]{interact}

\usepackage{epstopdf}% To incorporate .eps illustrations using PDFLaTeX, etc.
\usepackage[caption=false]{subfig}% Support for small, `sub' figures and tables
%\usepackage[nolists,tablesfirst]{endfloat}% To `separate' figures and tables from text if required
%\usepackage[doublespacing]{setspace}% To produce a `double spaced' document if required
%\setlength\parindent{24pt}% To increase paragraph indentation when line spacing is doubled

\usepackage[numbers,sort&compress]{natbib}% Citation support using natbib.sty
\bibpunct[, ]{[}{]}{,}{n}{,}{,}% Citation support using natbib.sty
\renewcommand\bibfont{\fontsize{10}{12}\selectfont}% Bibliography support using natbib.sty
\makeatletter% @ becomes a letter
\def\NAT@def@citea{\def\@citea{\NAT@separator}}% Suppress spaces between citations using natbib.sty
\makeatother% @ becomes a symbol again

\theoremstyle{plain}% Theorem-like structures provided by amsthm.sty
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{proposition}[theorem]{Proposition}

\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}

\theoremstyle{remark}
\newtheorem{remark}{Remark}
\newtheorem{notation}{Notation}


\usepackage{adjustbox, multirow}
\usepackage{amsmath}
\usepackage{amsfonts}

\providecommand{\norm}[1]{\lVert#1\rVert}
\DeclareMathOperator{\tr}{tr} 
\DeclareMathOperator{\littleo}{o}


\begin{document}

\articletype{ARTICLE TEMPLATE}% Specify the article type or omit as appropriate

\title{Supplementary Materials of the Paper Titled ``Comparing crossing hazard rate functions by joint modeling of survival and longitudinal data''}

\maketitle

\begin{abstract}
It is one of the important issues in survival analysis to compare two hazard rate functions to evaluate treatment effect. It is quite common that the two hazard rate functions cross each other at one or more unknown time points, representing temporal changes of the treatment effect. In certain applications, besides survival data, we also have related longitudinal data available regarding some time-dependent covariates. In such cases, a joint model that accommodates both types of data can allow us to infer the association between the survival and longitudinal data and to assess the treatment effect better. In this paper, we propose a modeling approach for comparing two crossing hazard rate functions by joint modeling survival and longitudinal data. Maximum likelihood estimation is used in estimating the parameters of the proposed joint model using the EM algorithm. Asymptotic properties of the maximum likelihood estimators are studied. To illustrate the virtues of the proposed method, we compare the performance of the proposed method with several existing methods in a simulation study. Our proposed method is also demonstrated using a real dataset obtained from an HIV clinical trial.
\end{abstract}

\begin{keywords}
censoring data; crossing hazard rates; longitudinal data; proportional hazards regression; survival analysis
\end{keywords}


\setcounter{section}{0}
\setcounter{equation}{0}
\def\theequation{S\arabic{section}.\arabic{equation}}
\def\thesection{S\arabic{section}}

\section{Proof of Lemma 1}
\setcounter{equation}{0}

Assume that $\lambda_{0}(u_{k}) \leq U$, for $k=1,\ldots, p$, where $U$ is a
finite positive number. From the assumption (C.2), the likelihood function 
$L$ is a continuous function of $\boldsymbol{\theta}$ on the compact set 
$Z \times \left[0,U\right]^{p}$. Then, $L$ can reach its maximum on the set 
by the extreme value theorem.

In order to show the existence of the MLE on $Z \times \left[0,\infty 
\right)^{p}$, it suffices to show that there exists a finite $U$ such that 
for each possible value of $\boldsymbol{\theta}_{U} \in \left\{Z \times 
\left[0,\infty \right)^{p}\right\}\setminus \left\{Z \times \left[0,U 
\right]^{p}\right\}$, there exists a $\boldsymbol{\theta} \in \left\{Z 
\times \left[0,U \right]^{p}\right\}$ such that $L(\boldsymbol{\theta}) 
> L(\boldsymbol{\theta}_{U})$. This can be done by using a proof of 
contradiction. Assume the existence of $\boldsymbol{\theta}_{U} \in 
\left\{Z \times \left[0,\infty \right)^{p}\right\}\setminus \left\{Z \times 
\left[0,U \right]^{p}\right\}$ that maximizes $L$ for each $U$. That is, 
$L(\boldsymbol{\theta}) < L(\boldsymbol{\theta}_{U})$ for $\boldsymbol{\theta} 
\in \left\{Z \times \left[0,U \right]^{p}\right\}$. Let
\[
   \min_{\boldsymbol{\zeta} \in Z, g\in \{0,1\}} {\exp\left\{\psi M \left( u \right)+\phi 
   \left(u-\gamma\right)g\right\}} = m > 0, 
\]
and
\[
   \max_{\boldsymbol{\zeta} \in Z, g\in \{0,1\}}  {\exp\left\{\psi M \left( u \right)+
   \phi \left(u-\gamma\right)g\right\}} = M > 0.
\]
The function $f_{3}^{'}\left(O_{i},\Delta_{i}| b_{i},\boldsymbol{\theta}_{t} 
\right)$ is then bounded above by
\[
   \left(M\prod_{k=1}^{p} \lambda_{0}\left(u_{k} \right)^{
   {I}{\left\{O_i=u_{k}\right\}}}\right)^{\Delta_{i}}  
   \exp{\left[-m\sum_{k=1}^{p}\lambda_{0}\left( u_{k} 
   \right){I}{\left\{O_i \geq u_{k}\right\}}\right]},
\]
and for convenience we denote this function by 
$f_{3}^{''}\left(O_{i},\Delta_{i}| \boldsymbol{b}_{i}, \boldsymbol{\theta}_{t} 
\right)$. Thus, the likelihood function $L$ is bounded above by
\[
	\prod_{i=1}^{n}\left[\int_{-\infty}^{\infty} 
                        \left\{ \prod_{j=1}^{n_i}f_{1}\left(Y(t_{ij}) | 
                        \boldsymbol{b}_{i}, \boldsymbol{\theta}_{y}\right) \right\} 
                        f_{2}\left(\boldsymbol{b}_{i} | \boldsymbol{\theta}_{b}\right)
                        f_{3}^{''}\left(O_{i},\Delta_{i}| \boldsymbol{b}_{i}, 
                        \boldsymbol{\theta}_{t}\right)d\boldsymbol{b}_{i}\right].
\]
From this result, we can show that $L(\boldsymbol{\theta}_{U})$ can be made 
as small as possible by increasing the value of $U$ as follows. If 
$\boldsymbol{\theta}_{U} \in \left\{Z \times \left[0,\infty \right)^{p}\right\}
\setminus \left\{Z \times \left[0,U \right]^{p}\right\}$, then there exists 
$1 \leq k \leq p$ such that $\lambda_{0}(u_{k})>U$. Thus, there exists at least 
one $1 \leq i \leq n$ such that $\sum_{k=1}^{p}\lambda_{0}\left( u_{k} \right){I}{\left\{O_i 
\geq u_{k}\right\}}>U$. Hence, $\sum_{k=1}^{p}\lambda_{0}\left( u_{k} \right)
{I}{\left\{O_i \geq u_{k}\right\}}$ goes to infinity as $U \rightarrow \infty$, 
and it follows that the upper bound of the likelihood $L(\boldsymbol{\theta}_{U})$ 
can be as close to $0$ as possible when $U$ increases. This is a contradiction
with the assumption that $\boldsymbol{\theta}_{U}$ maximizes the likelihood. 
Therefore, we can conclude that $\lambda_{0}(u_{k})$ are finite for all 
$k=1,\ldots, p$, and there is a finite $U$ such that the maximum likelihood is
within the compact set $Z \times \left[0,U\right]^{p}$ if it exists. Since the 
likelihood function is a continuous function of $\boldsymbol{\theta}$ on the 
compact set, the MLE $\hat{\boldsymbol{\theta}}$ exists.

From (8), the derivative of the likelihood function with respect to 
$\lambda_{0}(u)$ is
\[
   \frac{\sum_{i=1}^{n} \Delta_{i} I\left(O_{i}=u\right)}{\lambda_{0}(u)}
   - \sum_{i=1}^{n} E_{i}\left[\exp\left\{\psi M \left( u \right)+\phi 
   \left(u-\gamma\right)g_{i}\right\}\right]{I}{\left\{O_i \geq u\right\}}.
\]
After setting it to be 0, the maximum likelihood estimate of $\lambda_{0}(u)$ 
can be derived to be 
\[
   \hat{\lambda}_{0}(u)=\frac{\sum_{i=1}^{n} \Delta_{i} I\left(O_{i}=u\right)}
   {\sum_{i=1}^{n} E_{i}\left[\exp\left\{\psi M \left(u\right)+
   \phi \left(u-\gamma\right)g_{i}\right\}\right]{I}{\left\{O_i \geq u\right\}}}.
\]
The baseline hazard is then evaluated at each of the ordered failure times, 
denoted by $u_{1},u_{2},\ldots,u_{p}$. Then, the cumulative baseline hazard can be 
estimated by $\hat{\Lambda}(t)=\sum_{k=1}^{p} \hat{\lambda}_{0}
\left( u_{k} \right){I}{\left\{u_k \leq t\right\}}$, and it can be written as
\[
   \hat{\Lambda}(t)=\int_{0}^{t}\frac{dH_{n}(u)}
   {W_{n}(u;\hat{\boldsymbol{\theta}})},
\]
where $H_{n}(u)=\frac{1}{n}\sum_{i=1}^{n}\Delta_{i} I\left(O_{i}\leq u\right)$ 
and 
\[
   W_{n}(u;\hat{\boldsymbol{\theta}})=\frac{1}{n}\sum_{i=1}^{n} 
   E_{i}\big[\exp\{\hat{\psi} M \left(u\right)+\hat{\phi} 
   \left(u-\hat{\gamma}\right)g_{i}\}\big]{I}{\left\{O_i \geq u\right\}}.
\]


\section{Proof of Lemma 2}

Assume that there are $\boldsymbol{\theta}$ and $\boldsymbol{\theta}^{*}$ in 
$\Theta$ such that $L(\boldsymbol{\theta})=L(\boldsymbol{\theta}^{*})$ almost 
surely for all $n\geq 1$. We want to prove that
$\boldsymbol{\theta}=\boldsymbol{\theta}^{*}$. Since 
$L(\boldsymbol{\theta})=\prod_{i=1}^{n}L_{i}\left(\boldsymbol{\theta} | 
O_i,\Delta_i,\boldsymbol{Y_i},\boldsymbol{t_i}\right)$ and  
$L(\boldsymbol{\theta})=L(\boldsymbol{\theta}^{*})$ for all $n$, we have
$L_{i}\left(\boldsymbol{\theta} | O_i,\Delta_i,\boldsymbol{Y_i},
\boldsymbol{t_i}\right)=L_{i}\left(\boldsymbol{\theta}^{*} | 
O_i,\Delta_i,\boldsymbol{Y_i},\boldsymbol{t_i}\right)$, for all $1 \leq i \leq n$.
For simplicity, we omit the subscript $i$ in some expressions below. 
Now, $L_{i}\left(\boldsymbol{\theta} | O,\Delta,\boldsymbol{Y},\boldsymbol{t}\right)$ 
can be written as
\begin{equation}\label{eq:11}
\begin{split}
	&\int_{-\infty}^{\infty} \left(2\pi\sigma_{\varepsilon}^2\right)^{-\frac{N}{2}}
	\exp{\left\{-\frac{\left(\boldsymbol{Y}-\boldsymbol{X}'
        \boldsymbol{\beta}-\boldsymbol{Z}'\boldsymbol{b}\right)^{\intercal}
	\left(\boldsymbol{Y}-\boldsymbol{X}'\boldsymbol{\beta}-\boldsymbol{Z}'
        \boldsymbol{b}\right)}{2\sigma_{\varepsilon}^2}\right\}}\\
	&\times \left\{{\left(2\pi\right)^{q}|\Sigma_b|}\right\}^{-\frac{1}{2}}
	\exp{\left\{-\frac{1}{2}\boldsymbol{b}'\Sigma^{-1}_b\boldsymbol{b}\right\}}\\
	&\times  \lambda_{0}\left( O \right)^{\Delta}
	\left[\exp\left\{\psi\left(\boldsymbol{X}\left( O \right)'\boldsymbol{\beta}
	+\boldsymbol{Z}\left( O \right)'\boldsymbol{b}\right)
	+\phi\left(O - \gamma\right)\boldsymbol{g}\right\}\right]^{\Delta}\\
	&\times \exp{\left[-\int_{0}^{O}\lambda_{0}\left( u \right)
	\exp{\left\{\psi \left(\boldsymbol{X}\left( u \right)'\boldsymbol{\beta}
	+\boldsymbol{Z}\left( u \right)'\boldsymbol{b}\right)
	+\phi \left(u-\gamma\right)\boldsymbol{g}\right\}}du\right]}d\boldsymbol{b},
\end{split}
\end{equation}
where $N=\sum_{j=1}^n n_j$.
In the special case when $\Delta=0$ and $O=0$, the equation \eqref{eq:11} becomes
\[
	\int_{-\infty}^{\infty} \left(2\pi\sigma_{\varepsilon}^2\right)^{-\frac{N}{2}}
	\exp{\left\{-\frac{\left(\boldsymbol{Y}-\boldsymbol{X}'\beta-
        \boldsymbol{Z}'b\right)^{\intercal}\left(\boldsymbol{Y}-\boldsymbol{X}'
        \boldsymbol{\beta}-\boldsymbol{Z}'\boldsymbol{b}\right)}
        {2\sigma_{\varepsilon}^2}\right\}}
	\times \left\{{\left(2\pi\right)^{q}|\Sigma_b|}\right\}^{-\frac{1}{2}}
        \exp{\left\{-\frac{1}{2}\boldsymbol{b}'\Sigma^{-1}_b\boldsymbol{b}\right\}}
        d\boldsymbol{b},
\]
and it can be re-expressed as
\begin{equation}
\begin{split}
	&\left(2\pi\sigma_{\varepsilon}^2\right)^{-\frac{N}{2}}\frac{\left|
        \frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+
        \Sigma^{-1}_b\right|^{-\frac{1}{2}}}{\left|\Sigma_b\right|^{\frac{1}{2}}}  \\
	&\times \exp{\left\{-\frac{\left(\boldsymbol{Y}-\boldsymbol{X}'
        \boldsymbol{\beta}\right)'\left(\boldsymbol{Y}-\boldsymbol{X}'
        \boldsymbol{\beta} \right)}{2\sigma_{\varepsilon}^2}+ 
        \frac{1}{2\sigma_{\varepsilon}^4}\left(\boldsymbol{Y}-\boldsymbol{X}'
        \boldsymbol{\beta}\right)'\boldsymbol{Z}'\left(\frac{\boldsymbol{Z}
        \boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+\Sigma^{-1}_b\right)^{-1}
        \boldsymbol{Z}\left(\boldsymbol{Y}-\boldsymbol{X}'
        \boldsymbol{\beta}\right)\right\}}\\
	&\times \int_{-\infty}^{\infty} \left\{{\left(2\pi\right)^{q}
        \left|\frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+
        \Sigma^{-1}_b\right|^{-1}}\right\}^{-\frac{1}{2}}
	\times \exp\left\{-\frac{1}{2}\left[\boldsymbol{b}-
        \left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+
        \Sigma^{-1}_b\right)^{-1}\frac{\boldsymbol{Z}\left(\boldsymbol{Y}-
        \boldsymbol{X}'\boldsymbol{\beta}\right)}{\sigma_{\varepsilon}^2}\right]'\right. \\
 	&\left.\left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+
        \Sigma^{-1}_b\right) \left[\boldsymbol{b}-
        \left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+
        \Sigma^{-1}_b\right)^{-1}\frac{\boldsymbol{Z}\left(\boldsymbol{Y}-
        \boldsymbol{X}'\boldsymbol{\beta}\right)}{\sigma_{\varepsilon}^2}\right]\right\}
        d\boldsymbol{b}. 
        \nonumber
\end{split}
\end{equation}
Note that $\boldsymbol{b}$ is a random vector with a multivariate normal 
distribution with mean $\left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}
{\sigma_{\varepsilon}^2}+\Sigma^{-1}_b\right)^{-1}\frac{\boldsymbol{Z}
\left(\boldsymbol{Y}-\boldsymbol{X}'\boldsymbol{\beta}\right)}
{\sigma_{\varepsilon}^2}$ and covariance matrix $\left(\frac{\boldsymbol{Z}
\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+\Sigma^{-1}_b\right)^{-1}$. Thus, the 
equation \eqref{eq:11} becomes
\begin{equation}
\begin{split}
	&\left(2\pi\sigma_{\varepsilon}^2\right)^{-\frac{N}{2}}
	\frac{\left|\frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}
	+\Sigma^{-1}_b\right|^{-\frac{1}{2}}}{\left|\Sigma_b\right|^{\frac{1}{2}}}\\
	&\times \exp{\left\{-\frac{\left(\boldsymbol{Y}-\boldsymbol{X}'
        \boldsymbol{\beta}\right)'
	\left(\boldsymbol{Y}-\boldsymbol{X}'\boldsymbol{\beta} \right)}
	{2\sigma_{\varepsilon}^2}
	+ \frac{1}{2\sigma_{\varepsilon}^4}\left(\boldsymbol{Y}-\boldsymbol{X}'
        \boldsymbol{\beta}\right)'
	\boldsymbol{Z}'\left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}
	+\Sigma^{-1}_b\right)^{-1}\boldsymbol{Z}
	\left(\boldsymbol{Y}-\boldsymbol{X}'\boldsymbol{\beta}\right)\right\}}, \nonumber
\end{split}
\end{equation}
and this can be further written as
\begin{equation}
\begin{split}
	&\left(2\pi\sigma_{\varepsilon}^2\right)^{-\frac{N}{2}}
        \frac{\left|\frac{\boldsymbol{Z}\boldsymbol{Z}'}
        {\sigma_{\varepsilon}^2}+\Sigma^{-1}_b\right|^{-\frac{1}{2}}}
        {\left|\Sigma_b\right|^{\frac{1}{2}}}
	\times \exp \left\{-\frac{1}{2}\boldsymbol{Y}'\left[-\frac{1}
        {\sigma_{\varepsilon}^2}I_{q \times q} + \frac{1}{\sigma_{\varepsilon}^4}\boldsymbol{Z}'
        \left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+
        \Sigma^{-1}_b\right)^{-1}\boldsymbol{Z}\right]\boldsymbol{Y} \right.\\
	& \left.+ \boldsymbol{\beta}'\boldsymbol{X}\left[-\frac{1}
        {\sigma_{\varepsilon}^2}I_{q \times q} + \frac{1}{\sigma_{\varepsilon}^4}\boldsymbol{Z}'
        \left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+
        \Sigma^{-1}_b\right)^{-1}\boldsymbol{Z}\right]\boldsymbol{Y} \right.\\
	& \left.-\frac{1}{2} \boldsymbol{\beta}'\boldsymbol{X}\left[-\frac{1}
        {\sigma_{\varepsilon}^2}I_{q \times q} + \frac{1}{\sigma_{\varepsilon}^4}\boldsymbol{Z}'
        \left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+
        \Sigma^{-1}_b\right)^{-1}\boldsymbol{Z}\right]\boldsymbol{X}'
        \boldsymbol{\beta}\right\},\nonumber
\end{split}
\end{equation}
where $I_{q \times q}$ is the identity matrix with the same dimension as 
$\boldsymbol{Z}'\boldsymbol{Z}$.
Hence, under the conditions that $\Delta=0$, $O=0$, and 
$L_{i}\left(\boldsymbol{\theta} | O_i,\Delta_i,\boldsymbol{Y_i},
\boldsymbol{t_i}\right)=L_{i}\left(\boldsymbol{\theta}^{*} | 
O_i,\Delta_i,\boldsymbol{Y_i},\boldsymbol{t_i}\right)$ almost surely, we have
\begin{equation}
\begin{split}
	&\left(2\pi\sigma_{\varepsilon}^2\right)^{-\frac{N}{2}}
        \frac{\left|\frac{\boldsymbol{Z}\boldsymbol{Z}'}
        {\sigma_{\varepsilon}^2}+\Sigma^{-1}_b\right|^{-\frac{1}{2}}}
        {\left|\Sigma_b\right|^{\frac{1}{2}}}
	\times \exp \left\{-\frac{1}{2}\boldsymbol{Y}'
        \left[-\frac{1}{\sigma_{\varepsilon}^2}I_{q \times q} + \frac{1}{\sigma_{\varepsilon}^4}
        \boldsymbol{Z}'\left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}
        {\sigma_{\varepsilon}^2}+\Sigma^{-1}_b\right)^{-1}\boldsymbol{Z}\right]
        \boldsymbol{Y} \right.\\
	& \left.+ \boldsymbol{\beta}'\boldsymbol{X}\left[-\frac{1}
        {\sigma_{\varepsilon}^2}I_{q \times q} + \frac{1}{\sigma_{\varepsilon}^4}\boldsymbol{Z}'
        \left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+
        \Sigma^{-1}_b\right)^{-1}\boldsymbol{Z}\right]\boldsymbol{Y} \right.\\
	& \left.-\frac{1}{2} \boldsymbol{\beta}'\boldsymbol{X}\left[-\frac{1}
        {\sigma_{\varepsilon}^2}I_{q \times q} + \frac{1}{\sigma_{\varepsilon}^4}\boldsymbol{Z}'
        \left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+
        \Sigma^{-1}_b\right)^{-1}\boldsymbol{Z}\right]\boldsymbol{X}'
        \boldsymbol{\beta}\right\}\\
	&=\left(2\pi{\sigma_{\varepsilon}^{*}}^{2}\right)^{-\frac{N}{2}}
	\frac{\left|\frac{\boldsymbol{Z}\boldsymbol{Z}'}
        {{\sigma_{\varepsilon}^{*}}^{2}}+{\Sigma^{*}_b}^{-1}\right|^{-\frac{1}{2}}}
	{\left|{\Sigma^{*}_b}\right|^{\frac{1}{2}}}
	\times \exp \left\{-\frac{1}{2}\boldsymbol{Y}'\left[-\frac{1}
        {{\sigma_{\varepsilon}^{*}}^{2}}I_{q \times q} + \frac{1}{{\sigma_{\varepsilon}^{*}}^{4}}
        \boldsymbol{Z}'\left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}
        {{\sigma_{\varepsilon}^{*}}^{2}}+{\Sigma^{*}_b}^{-1}\right)^{-1}
        \boldsymbol{Z}\right]\boldsymbol{Y} \right.\\
	& \left.+ \boldsymbol{\beta^{*}}'\boldsymbol{X}\left[-\frac{1}
        {{\sigma_{\varepsilon}^{*}}^{2}}I_{q \times q} + \frac{1}{{\sigma_{\varepsilon}^{*}}^{4}}
        \boldsymbol{Z}'\left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}
        {{\sigma_{\varepsilon}^{*}}^{2}}+{\Sigma^{*}_b}^{-1}\right)^{-1}
        \boldsymbol{Z}\right]\boldsymbol{Y} \right.\\
	& \left.-\frac{1}{2} \boldsymbol{\beta^{*}}'\boldsymbol{X}
        \left[-\frac{1}{{\sigma_{\varepsilon}^{*}}^{2}}I_{q \times q} + \frac{1}
        {{\sigma_{\varepsilon}^{*}}^{4}}\boldsymbol{Z}'
        \left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{{\sigma_{\varepsilon}^{*}}^{2}}+
        {\Sigma^{*}_b}^{-1}\right)^{-1}\boldsymbol{Z}\right]
        \boldsymbol{X}'\boldsymbol{\beta^{*}}\right\}.\nonumber
\end{split}
\end{equation}
Then, we get
\begin{equation}
\begin{split}
	\frac{1}{\sigma_{\varepsilon}^N}
	\frac{\left|\frac{\boldsymbol{Z}\boldsymbol{Z}'}
	{\sigma_{\varepsilon}^2}+\Sigma_b^{-1}\right|^{-\frac{1}{2}}}
	{\left|\Sigma_b\right|^{\frac{1}{2}}} 
	& =\frac{1}{{\sigma_{\varepsilon}^{*}}^N}
	\frac{\left|\frac{\boldsymbol{Z}\boldsymbol{Z}'}
	{{\sigma_{\varepsilon}^{*}}^2}+{\Sigma^{*}_b}^{-1}\right|^{-\frac{1}{2}}}
	{\left|{\Sigma^{*}_b}\right|^{\frac{1}{2}}} \\
	\boldsymbol{X}'\boldsymbol{\beta}  
	&= \boldsymbol{X}'\boldsymbol{\beta^{*}}\nonumber
\end{split}
\end{equation}
and
\begin{equation}\label{eq:12}
	- \frac{1}{\sigma_{\varepsilon}^2}I_{q \times q} 
	+ \frac{1}{\sigma_{\varepsilon}^4}\boldsymbol{Z}'
        \left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}
	+\Sigma_b^{-1}\right)^{-1}\boldsymbol{Z}
 	=- \frac{1}{{\sigma_{\varepsilon}^{*}}^2}I_{q \times q} 
	+ \frac{1}{{\sigma_{\varepsilon}^{*}}^4}\boldsymbol{Z}'
        \left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{{\sigma_{\varepsilon}^{*}}^2}
	+{\Sigma^{*}_b}^{-1}\right)^{-1}\boldsymbol{Z}.
\end{equation}

After some algebraic manipulations and by the assumption (C.4), the equation 
\eqref{eq:12} becomes
\[
	- \frac{1}{\sigma_{\varepsilon}^2}I_{q \times q} 
	+ \frac{1}{\sigma_{\varepsilon}^4}\left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}
        {\sigma_{\varepsilon}^2}
	+\Sigma_b^{-1}\right)^{-1}\boldsymbol{Z}\boldsymbol{Z}'
	=- \frac{1}{{\sigma_{\varepsilon}^{*}}^2}I_{q \times q} + \frac{1}{{\sigma_{\varepsilon}^{*}}^4}
        \left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}{{\sigma_{\varepsilon}^{*}}^2}
	+{\Sigma^{*}_b}^{-1}\right)^{-1}\boldsymbol{Z}\boldsymbol{Z}',
\]
and the left-hand-side of this equation is 
\begin{equation}
\begin{split}
	- \frac{1}{\sigma_{\varepsilon}^2}I_{q \times q} 
	+ \frac{1}{\sigma_{\varepsilon}^4}\left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}
        {\sigma_{\varepsilon}^2}
	+\Sigma_b^{-1}\right)^{-1}\boldsymbol{Z}\boldsymbol{Z}'
	&=\frac{1}{\sigma_{\varepsilon}^4}\left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}
        {\sigma_{\varepsilon}^2}
	+\Sigma_b^{-1}\right)^{-1} \left[-\sigma_{\varepsilon}^2\left(
        \frac{\boldsymbol{Z}\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}
	+\Sigma_b^{-1}\right)+\boldsymbol{Z}\boldsymbol{Z}'\right]\\
	&= \frac{1}{\sigma_{\varepsilon}^4}\left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}
        {\sigma_{\varepsilon}^2}
	+\Sigma_b^{-1}\right)^{-1} \left[-\sigma_{\varepsilon}^2 \Sigma_b^{-1}\right]\\
	&= - \left(\boldsymbol{Z}\boldsymbol{Z}'+\sigma_{\varepsilon}^2
        \Sigma_b^{-1}\right)^{-1} \Sigma_b^{-1}\\
	&= - \left(\Sigma_b\boldsymbol{Z}\boldsymbol{Z}'
	+\sigma_{\varepsilon}^2I_{q \times q}\right)^{-1}.\nonumber
\end{split}
\end{equation}
Similarly, the right-hand-side is $- \left(\Sigma_b^{*}\boldsymbol{Z}\boldsymbol{Z}'
+\sigma^{*2}_{\varepsilon}I_{q \times q}\right)^{-1}$. Thus, we obtain
\[
   \Sigma_b \boldsymbol{Z}\boldsymbol{Z}' + \sigma_{\varepsilon}^2 I_{q \times q} = 
   {\Sigma_b^{*}} \boldsymbol{Z}\boldsymbol{Z}' + {\sigma^{*}_{\varepsilon}}^2 I_{q \times q} . 
\]
By combining all these results and the assumption (C.4), we obtain 
$\boldsymbol{\beta}=\boldsymbol{\beta^{*}}$, $\sigma_{\varepsilon}^2=
{\sigma^{*}_{\varepsilon}}^2$, and $\Sigma_b=\Sigma^{*}_b$.

When $\Delta=0$, the equation \eqref{eq:11} becomes
\begin{equation}\label{eq:13}
\begin{split}
	&\int_{-\infty}^{\infty} \left(2\pi\sigma_{\varepsilon}^2\right)^{-\frac{N}{2}}
	\exp{\left\{-\frac{\left(\boldsymbol{Y}-\boldsymbol{X}'
        \boldsymbol{\beta}-\boldsymbol{Z}'\boldsymbol{b}\right)'
	\left(\boldsymbol{Y}-\boldsymbol{X}'\boldsymbol{\beta}-
        \boldsymbol{Z}'\boldsymbol{b}\right)}
	{2\sigma_{\varepsilon}^2}\right\}}\\
	&\times \left\{{\left(2\pi\right)^{q}|\Sigma_b|}\right\}^{-\frac{1}{2}}
	\exp{\left\{-\frac{1}{2}\boldsymbol{b}'\Sigma_b^{-1} \boldsymbol{b}\right\}}\\
	&\times \exp\left[-\int_{0}^{O}\lambda_{0}\left( u \right)
	\exp\left\{\psi \left(\boldsymbol{X}\left( u \right)'
        \beta+\boldsymbol{Z}\left( u \right)'\boldsymbol{b}\right)
	+\phi \left(u-\gamma\right)\boldsymbol{g}\right\}du\right]d\boldsymbol{b}.
\end{split}
\end{equation}
As before, $\boldsymbol{b}$ is a random vector with the multivariate normal 
distribution with mean $\left(\frac{\boldsymbol{Z}\boldsymbol{Z}'}
{\sigma_{\varepsilon}^2}+\Sigma_b^{-1}\right)^{-1}\frac{\boldsymbol{Z}
\left(\boldsymbol{Y}-\boldsymbol{X}'\boldsymbol{\beta}\right)}
{\sigma_{\varepsilon}^2}$ and covariance matrix $\left(\frac{\boldsymbol{Z}
\boldsymbol{Z}'}{\sigma_{\varepsilon}^2}+\Sigma_b^{-1}\right)^{-1}$. After applying
this distribution to the equation \eqref{eq:13}, the equation can be regarded as 
the expectation of 
\[
	\exp{\left[-\int_{0}^{O}\lambda_{0}\left( u \right)
	\exp{\left\{\psi \left(\boldsymbol{X}\left( u \right)'\boldsymbol{\beta}
	+\boldsymbol{Z}\left( u \right)'\boldsymbol{b}\right)
	+\phi \left(u-\gamma\right)\boldsymbol{g}\right\}}du\right]}.
\]
For fixed $\boldsymbol{X}$ and $\boldsymbol{Z}$, $\boldsymbol{b}$ is the 
complete statistics for $\boldsymbol{Z}\boldsymbol{Y}$ by considering 
$\boldsymbol{Z}\boldsymbol{Y}$ as a parameter. Hence, under $\Delta=0$, 
$L_{i}\left(\boldsymbol{\theta} | O_i,\Delta_i,\boldsymbol{Y_i},
\boldsymbol{t_i}\right)=L_{i}\left(\boldsymbol{\theta}^{*} | 
O_i,\Delta_i,\boldsymbol{Y_i},\boldsymbol{t_i}\right)$ almost surely implies 
that
\begin{equation}
\begin{split}
	&\exp{\left[-\int_{0}^{O}\lambda_{0}\left( u \right)
	\exp{\left\{\psi \left(\boldsymbol{X}\left( u \right)'\boldsymbol{\beta}
	+\boldsymbol{Z}\left( u \right)'\boldsymbol{b}\right)
	+\phi \left(u-\gamma\right)\boldsymbol{g}\right\}}du\right]}\\
	&=\exp{\left[-\int_{0}^{O}\lambda_{0}^{*}\left( u \right)
	\exp{\left\{\psi^{*} \left(\boldsymbol{X}\left( u \right)'\boldsymbol{\beta^{*}}
	+\boldsymbol{Z}\left( u \right)'\boldsymbol{b}\right)
	+\phi^{*} \left(u-\gamma^{*}\right)\boldsymbol{g}\right\}}du\right]},\nonumber
\end{split}
\end{equation}
and this implies that
\begin{equation}
\begin{split}
	&\lambda_{0}\left( u \right)
	\exp{\left\{\psi \left(\boldsymbol{X}\left( u \right)'\boldsymbol{\beta}
	+\boldsymbol{Z}\left( u \right)'\boldsymbol{b}\right)
	+\phi \left(u-\gamma\right)\boldsymbol{g}\right\}}\\
	&=\lambda_{0}^{*}\left( u \right)
	\exp{\left\{\psi^{*} \left(\boldsymbol{X}\left( u \right)'\boldsymbol{\beta^{*}}
	+\boldsymbol{Z}\left( u \right)'\boldsymbol{b}\right)
	+\phi^{*} \left(u-\gamma^{*}\right)\boldsymbol{g}\right\}}.\nonumber
\end{split}
\end{equation}
By the assumptions (C.2) and (C.5), we have that $\psi={\psi^{*}}$, 
$\boldsymbol{\beta}=\boldsymbol{\beta^{*}}$, $\phi={\phi^{*}}$, 
$\gamma={\gamma^{*}}$, and $\Lambda_0={\Lambda^{*}_0}$. Hence, we can conclude 
that the MLEs of the parameters are unique under the assumptions (C.1) - (C.5).


\section{Proof of Theorem 1}

The proof has two parts. First, we will show that 
$\hat{\Lambda}_0(t)$ is bounded almost surely as $n \rightarrow \infty$. 
Then, we will show that any 
convergent subsequence of $\hat{\boldsymbol{\theta}}=(\hat{\boldsymbol{\zeta}}, 
\hat{\Lambda}_0)$ must converge to the true value $\boldsymbol{\theta}=
(\boldsymbol{\zeta}, \Lambda_{0})$, using the identifiability result. Next, we
justify these results in details.

% The first part %
First, for all $t \in \left[0,\tau\right]$, we have $\hat{\Lambda}_0(t)\geq0$ and 
\[
	\hat{\Lambda}_0(\tau) 
	\leq \frac{\sum_{i=1}^{n} \Delta_{i} I\left(O_{i}\leq \tau\right)}
	{m\sum_{i=1}^{n}I\left(O_{i}\geq \tau\right)}.
\]
Note that $\frac{1}{n}\sum_{i=1}^{n} \Delta_{i} I\left(O_{i}\leq \tau\right) 
\leq\frac{1}{n}\sum_{i=1}^{n} I\left(O_{i}\leq \tau\right)$, 
$\frac{1}{n}\sum_{i=1}^{n} I\left(O_{i}\leq \tau\right)$ converges to 
$E_{0}\left(O\leq \tau\right)$ as $n \rightarrow \infty$, and that 
$\frac{1}{n}\sum_{i=1}^{n} I\left(O_{i}\geq \tau\right)$ converges to 
$E_{0}\left(O\geq \tau\right)$ as $n \rightarrow \infty$. Thus, 
$\hat{\Lambda}_0(\tau)$ stays bounded as $n \rightarrow \infty$. 
Since the cumulative hazard rate function is a nondecreasing function 
and the bound does not depend on $t$, we conclude that $\hat{\Lambda}_{0}$ 
is uniformly bounded almost surely.

% The second part %
For the second part of the proof, we will repeatedly use the following result. 
The class $\boldsymbol{BV}_{\left[0,\tau\right]}$ of all functions $f:\left[0,\tau\right] 
\rightarrow \mathbb{R}$ that are uniformly bounded by a constant $M$ and of 
variation bounded by $M$ is Donsker. Also, every uniformly bounded Donsker 
class is Glivenko-Cantelli; that is, $\sup_{f \in \mathbb{F}} |\mathbb{E}_{n}f - E_{0}f| 
\rightarrow 0$ almost surely. More details of these results can be found in 
\cite{van1996weak}.

By the Bolzanno-Weierstrass theorem, each bounded sequence 
$(\hat{\boldsymbol{\zeta}})_{n \in \mathbb{N}}$ in $\mathbb{R}^r$ has a 
convergent subsequence, which converges to some $\boldsymbol{\zeta^{*}}$. 
By the Glivenko-Cantelli theorem, $\left(H_{n}(u)\right)_{n \in \mathbb{N}}$ converges 
uniformly to $H(u)$ on $\left[0,\tau\right]$. Since $\hat{\Lambda}_0$ is 
bounded on $\left[0,\tau\right]$, the functions 
$u \rightarrow W_{n}(u;\hat{\boldsymbol{\theta}})$ are uniformly bounded 
away from both $0$ and $\infty$. By the Helly-Bray lemma, we can choose a 
subsequence of $W_{n}(u;\hat{\boldsymbol{\theta}})$ such that the subsequence 
converges to some monotone function $W^{*}(u)$ for every $u$. Then, by the 
dominated convergence theorem, $\hat{\Lambda}_0$ converges uniformly to 
$\Lambda_0^{*}=\int_{0}^{t}\frac{dH(u)}{W^{*}(u)}$ on $\left[0,\tau\right]$. More 
specifically, 
\[
	\hat{\Lambda}_0(t)
	=\int_{0}^{t}\frac{dH_{n}(u)}{W_{n}(u;\hat{\boldsymbol{\theta}})}
	=\int_{0}^{t}\frac{dH(u)}{W_{n}(u;\hat{\boldsymbol{\theta}})} + o(1) 
	\rightarrow  \Lambda_0^{*}(t)=\int_{0}^{t}\frac{dH(u)}{W^{*}(u)}
\]
uniformly for $t \in \left[0,\tau\right]$. Hence, 
$\norm{\hat{\Lambda}_0-\Lambda_0^{*}}_\infty$ converges almost surely to 0.
Let $\bar{\Lambda}_0(t)=\int_{0}^{t}\frac{dH_{n}(u)}{W_{n}(u;\boldsymbol{\theta})}$ 
be the intermediate quantity between $\hat{\Lambda}_0(t)=\int_{0}^{t}
\frac{dH_{n}(u)}{W_{n}(u;\hat{\boldsymbol{\theta}})}$ and 
$\Lambda_{0}\left(t\right)=\int_{0}^{t}\frac{dH(u)}{W(u;\boldsymbol{\theta})}$ 
where $W(u;\boldsymbol{\theta})$ and $H(u)$ are the expectations of 
$W_{n}(u;\boldsymbol{\theta})$ and $H_{n}(u)$, evaluated at $\boldsymbol{\theta}$. 
That is, $H(u)=E_{0}\left[\Delta I\left(O\leq u\right)\right]$ and 
$$W(u;\boldsymbol{\theta})=E_{0}\left[\exp\left\{\psi M \left(u\right)+\phi 
\left(u-\gamma\right)g\right\}{1}{\left\{O \geq u\right\}}\right].$$ Next, we
prove that $\bar{\Lambda}_{0}\left(t\right)$ converges uniformly to 
$\Lambda_{0}\left(t\right)$ on $\left[0,\tau\right]$. 
By the Glivenko-Cantelli theorem, $\left(H_{n}(u)\right)_{n \in \mathbb{N}}$ converges 
uniformly to $H(u)$ on $\left[0,\tau\right]$. Since $\Lambda_{0}$ is bounded on 
$\left[0,\tau\right]$ and $\boldsymbol{\zeta}$ is in a compact set, the functions 
$u \rightarrow W_{n}(u;\boldsymbol{\theta})$ are contained in $\boldsymbol{BV}_{\left[0,\tau\right]}$ and uniformly bounded away from $0$. Thus, by the Glivenko-Cantelli theorem, 
$\left(W_{n}(u;\boldsymbol{\theta})\right)_{n \in \mathbb{N}}$ converges uniformly  
on $\left[0,\tau\right]$ to $W(u;\boldsymbol{\theta})$, which is bounded away 
from $0$. Hence, by the dominated convergence theorem, 
$\norm{\bar{\Lambda}_{0}-\Lambda_0}_\infty$ converges almost surely to $0$.

% The third part %
Finally, we will show that any convergent subsequence of 
$\hat{\boldsymbol{\theta}}$ must converge to the true value 
$\boldsymbol{\theta}$, using the identifiability result. From the second part of 
the proof, we have 
\[
	\log{L(\hat{\boldsymbol{\zeta}}, \hat{\Lambda}_0)}
	-\log{L(\boldsymbol{\zeta}, \bar{\Lambda}_0)} 
	\rightarrow \log{L(\boldsymbol{\zeta^{*}},\Lambda_0^{*})}
	-\log{L(\boldsymbol{\zeta}, \Lambda_{0})}.
	\qquad \mbox{a.s.}
\]
By the dominated convergence theorem,
\[
	E_{0}\left[\log{L(\hat{\boldsymbol{\zeta}}, \hat{\Lambda}_0)}
	-\log{L(\boldsymbol{\zeta}, \bar{\Lambda}_0)}\right] 
	\rightarrow E_{0}\left[\log{L(\boldsymbol{\zeta^{*}},\Lambda_0^{*})}
	-\log{L(\boldsymbol{\zeta}, \Lambda_{0})}\right]. 
	\qquad \mbox{a.s.}	
\]
Because $\hat{\boldsymbol{\theta}}=(\hat{\boldsymbol{\zeta}}, 
\hat{\Lambda}_0)$ maximizes the log-likelihood,
\[
	\frac{1}{n}\sum_{i=1}^{n} 
	\left[\log{L_{i}(\hat{\boldsymbol{\zeta}}, \hat{\Lambda}_0)}
	-\log{L_{i}(\boldsymbol{\zeta}, \bar{\Lambda}_0)}\right] \geq 0 .
\]
It follows that as $n \rightarrow \infty$,
\[
	E_{0} \left[\log{L(\hat{\boldsymbol{\zeta}}, \hat{\Lambda}_0)}
	-\log{L(\boldsymbol{\zeta}, \bar{\Lambda}_0)}\right] \geq 0.
	\qquad \mbox{a.s.}	
\]
Hence, $E_{0}\left[\log{L(\boldsymbol{\zeta^{*}},\Lambda_0^{*})}-
\log{L(\boldsymbol{\zeta}, \Lambda_{0})}\right]  \geq 0$. However, the left-hand 
side of the inequality is the negative Kullback-Leibler information; that 
is, $E_{0}\left[\log{L(\boldsymbol{\zeta^{*}},\Lambda_0^{*})}-
\log{L(\boldsymbol{\zeta}, \Lambda_{0})}\right]  \leq 0$. Thus,  
$E_{0}\left[\log{L(\boldsymbol{\zeta^{*}},\Lambda_0^{*})}-
\log{L(\boldsymbol{\zeta}, \Lambda_{0})}\right]  = 0$, and this implies that 
$L(\boldsymbol{\zeta^{*}},\Lambda_0^{*})=L(\boldsymbol{\zeta}, \Lambda_{0})$ 
almost everywhere. Hence, $\boldsymbol{\zeta^{*}}=\boldsymbol{\zeta}$ and 
$\Lambda_0^{*}=\Lambda_{0}$ by the identifiability result.


\section{Proof of Theorem 2}

The proof of the asymptotic normality of the parameter estimators of our proposed
joint models is based on the Theorem 3.3.1 of \cite{van1996weak}. 
We will verify that the conditions of that theorem are all satisfied in the current 
setting. Then, the asymptotic normality of our estimators follows directly.

The empirical and expected scores are used in the proof. First, we show that they 
are given by
\begin{eqnarray}\label{eq:14}
	S_{n}\left(\boldsymbol{\theta}\right)\left(\boldsymbol{h_{1}},h_{2}\right)
	&=&\mathbb{E}_{n}\left\{l_{1}\left(\boldsymbol{\theta}\right)'\boldsymbol{h_{1}}
	+ l_{2}\left(\boldsymbol{\theta}\right)[h_{2}]\right\} \nonumber\\
	S\left(\boldsymbol{\theta}\right)\left(\boldsymbol{h_{1}},h_{2}\right)
	&=&E_{0}\left\{l_{1}\left(\boldsymbol{\theta}\right)'\boldsymbol{h_{1}} 
	+ l_{2}\left(\boldsymbol{\theta}\right)[h_{2}]\right\},
\end{eqnarray}
where $l_{1}\left(\boldsymbol{\theta}\right)'\boldsymbol{h_{1}}$ is the score 
function for $\boldsymbol{\zeta}$ in the direction of $\boldsymbol{h_{1}}$ and 
$ l_{2}\left(\boldsymbol{\theta}\right)[h_{2}]$ is the score function for 
$\Lambda_0$ in the direction of $h_{2}$. To obtain the score functions, let us 
consider the one-dimensional submodels through the estimators, differentiate 
them with respect to the parameters and evaluate the differentiations at the
estimators, instead of directly differentiating the log-likelihood function with 
respect to $\boldsymbol{\zeta}$ and the jump sizes of $\Lambda_0$. More specifically, 
let $\boldsymbol{h}$ be in the space of $H=\{\boldsymbol{h}=
\left(\boldsymbol{h_{1}},h_{2}\right)|\norm{ \boldsymbol{h_{1}}} \leq 1, 
\norm{h_2}_{v} \leq 1\}$ where $\norm{ \cdot }$ 
is the Euclidean norm and $\norm{h_{2}}_{v}$ is the absolute value of $h_{2}(0)$ plus 
the total variation of $h_{2}$ on $\left[0,\tau\right]$. We then define 
$\boldsymbol{\zeta}_{t}=\boldsymbol{\zeta}+t \boldsymbol{h_{1}}$, 
$\Lambda_{t}(\cdot)=\int_{0}^{\cdot}\left(1+t h_{2}(u)\right)d\Lambda_0(u)$, and 
$\theta_{t}=\left(\boldsymbol{\zeta}_{t},\Lambda_{t}\right)$, for $\boldsymbol{h} \in H$. 
Then, we can obtain the score $S_{n}\left(\boldsymbol{\theta} \right)
\left(\boldsymbol{h}\right)$ after differentiating with respect to $t$ and 
evaluating at $t=0$.
For $\boldsymbol{h}=\left(\boldsymbol{h_{1}},h_{2}\right) \in H$, let 
$\left(h_{1}^{\boldsymbol{\beta}},h_{1}^{\sigma_{\varepsilon}^2},h_{1}^{\Sigma_b},h_{1}^{\psi},
h_{1}^{\phi},h_{1}^{\gamma}\right)$ be the corresponding components of $\boldsymbol{h_{1}}$ 
for $\left(\boldsymbol{\beta},\sigma_{\epsilon}^2,\Sigma_b,\psi,\phi,\gamma \right)$, 
respectively. Because the log-likelihood is maximized at 
$(\hat{\boldsymbol{\zeta}},\hat{\Lambda}_0)$, it is maximized at $t=0$ when 
evaluated on the submodel $(\boldsymbol{\zeta}_{t},\hat{\Lambda}_0)$.
Thus, we can obtain the score function for $\boldsymbol{\zeta}$ after 
differentiating the log-likelihood with respect to $t$ and evaluating at $t=0$. 
Similarly, we can obtain the score function for $\Lambda_0$ through the submodel 
$(\hat{\boldsymbol{\zeta}},\Lambda_{t})$.

For simplicity, let us denote the observed data by $D_{o}$. Then, the log-likelihood 
function $l(\boldsymbol{\theta}_{t})$ for a subject can be written as
\begin{equation}
\begin{split}
	&E\left[-\frac{N}{2}\log{\left\{2\pi(\sigma_{\varepsilon}^2+
        th_{1}^{\sigma_{\varepsilon}^2})\right\}}
	-\frac{\left(\boldsymbol{Y}-\boldsymbol{X}'(\boldsymbol{\beta}+
        th_{1}^{\boldsymbol{\beta}})-\boldsymbol{Z}'\boldsymbol{b}\right)'
	\left(\boldsymbol{Y}-\boldsymbol{X}'(\boldsymbol{\beta}+
        th_{1}^{\boldsymbol{\beta}})-\boldsymbol{Z}'\boldsymbol{b}\right)}
	{2(\sigma_{\varepsilon}^2+th_{1}^{\sigma_{\varepsilon}^2})}\right.\\
	&\left. -\frac{1}{2} \log{\left\{{\left(2\pi\right)^{q}\left|
        \Sigma_b + t h_{1}^{\Sigma_b}\right|}\right\}} 
	-\frac{1}{2}\boldsymbol{b}' \left(\Sigma_b + t h_{1}^{\Sigma_b}\right)^{-1} 
        \boldsymbol{b}
	+ \Delta \log{\left\{\left( 1 + t h_{2}(O) \right)d\Lambda_0(O)\right\}}\right.\\
	&\left. + \Delta \left\{\left(\psi + t h_{1}^\psi\right)
	\left(\boldsymbol{X}\left( O \right)'(\boldsymbol{\beta}+th_{1}^{\boldsymbol{\beta}})
	+\boldsymbol{Z}\left( O \right)'\boldsymbol{b}\right)
	+\left(\phi + t h_{1}^\phi\right)\left(O - \left(\gamma + 
        th_{1}^\gamma\right) \right)\boldsymbol{g}\right\}\right.\\
	&\left. - \int_{0}^{O} \exp\Big\{\left(\psi + t h_{1}^\psi\right) 
	\left(\boldsymbol{X}\left( u \right)'(\boldsymbol{\beta}+th_{1}^{\boldsymbol{\beta}})
	+\boldsymbol{Z}\left( u \right)'\boldsymbol{b}\right)\right.\\
	&\left. +\left(\phi + t h_{1}^\phi\right)\left(u-\left(\gamma + 
        th_{1}^\gamma\right)\right)\boldsymbol{g}\Big\} 
	\left( 1 + t h_{2}(u) \right)d\Lambda\left(u\right)\big|D_{o}\right].\nonumber
\end{split}
\end{equation}
Then, by differentiating it with respect to $t$ and evaluating the differentiation 
at $t=0$, we get
\begin{equation}
\begin{split}
	\frac{\partial l(\boldsymbol{\theta}_{t})}{\partial t}\bigg|_{t=0} 
	& =:
	l_{1}\left(\boldsymbol{\theta}\right)'\boldsymbol{h_{1}} 
	+ l_{2}\left(\boldsymbol{\theta}\right)[h_{2}]\\
	&=m_{1}\left(O | \boldsymbol{\theta}\right)'\boldsymbol{h_{1}} 
	- \int_{0}^{O}m_{2}\left(u | \boldsymbol{\theta}\right)'\boldsymbol{h_{1}} d\Lambda_0(u)
	+\Delta h_{2}(O) 
	- \int_{0}^{O}m_{3}\left(u | \boldsymbol{\theta}\right)h_{2}(u) d\Lambda_0(u),\nonumber
\end{split}
\end{equation}
where
\begin{equation}
\begin{split}
	m_{1}\left(O | \boldsymbol{\theta}\right)'\boldsymbol{h_{1}}
 	&=E\left[-\frac{Nh_{1}^{\sigma_{\varepsilon}^2}}{2\sigma_{\varepsilon}^2}
	+\frac{\left(\boldsymbol{Y}-\boldsymbol{X}'\boldsymbol{\beta}
	-\boldsymbol{Z}'\boldsymbol{b}\right)'\left(\boldsymbol{Y}-
        \boldsymbol{X}'\boldsymbol{\beta}-\boldsymbol{Z}'\boldsymbol{b}\right)
	h_{1}^{\sigma_{\varepsilon}^2}}{\sigma_{\varepsilon}^4}\right.\\
	&\left.+\frac{\boldsymbol{X}\left(\boldsymbol{Y}-
        \boldsymbol{X}'\boldsymbol{\beta}-\boldsymbol{Z}'\boldsymbol{b}\right)
	h_{1}^{\boldsymbol{\beta}}}{\sigma_{\varepsilon}^2}
	-\frac{\boldsymbol{b}'\Sigma_b^{-1}h_{1}^{\Sigma_b}\Sigma_b^{-1}\boldsymbol{b}}{2} 
        - \tr{\left(\Sigma_b^{-1}h_{1}^{\Sigma_b}\right)}\right.\\
	&+ \Delta \left\{h_{1}^\psi\left(\boldsymbol{X}\left( O \right)'\boldsymbol{\beta}
	+\boldsymbol{Z}\left( O \right)'\boldsymbol{b}\right)
	+ \psi\left(\boldsymbol{X}\left( O \right)'h_{1}^{\boldsymbol{\beta}}
	+\boldsymbol{Z}\left( O \right)'\boldsymbol{b}\right) \right. \\
	& \left.\left. + h_{1}^\phi \left(O - \gamma \right)\boldsymbol{g} 
	+ \phi \left(O - h_{1}^\gamma \right)\boldsymbol{g}\right\}\big|D_{o}\right],\nonumber
\end{split}
\end{equation}
\begin{equation}
\begin{split}
	m_{2}\left(u | \boldsymbol{\theta}\right)'\boldsymbol{h_{1}}
	&=E\left[\left\{h_{1}^\psi\left(\boldsymbol{X}\left( u \right)'\boldsymbol{\beta}
	+\boldsymbol{Z}\left( u \right)'\boldsymbol{b}\right)
	+ \psi\left(\boldsymbol{X}\left( u \right)'h_{1}^{\boldsymbol{\beta}}
	+\boldsymbol{Z}\left( u \right)'\boldsymbol{b}\right)+ \right.\right.\\
        & \left. h_{1}^\phi \left(u - \gamma \right)\boldsymbol{g} 
	+ \phi \left(u - h_{1}^\gamma \right)\boldsymbol{g}\right\} \\
	&\left. \times \exp\left\{\psi\left(\boldsymbol{X}\left( u \right)'
        \boldsymbol{\beta}+\boldsymbol{Z}\left( u \right)'\boldsymbol{b}\right)+
        \phi\left(u-\gamma\right)\boldsymbol{g}\right\}\big|D_{o}\right],\nonumber
\end{split}
\end{equation}
and
\[
	m_{3}\left(u | \boldsymbol{\theta}\right)
	=E\left[\exp\left\{\psi\left(\boldsymbol{X}\left( u \right)'\boldsymbol{\beta}
	+\boldsymbol{Z}\left( u \right)'\boldsymbol{b}\right)
	+\phi\left(u-\gamma\right)\boldsymbol{g}\right\}\big|D_{o}\right].
\]
Hence, the empirical and expected scores are given by \eqref{eq:14} .
Let $\boldsymbol{\theta}(\boldsymbol{h})=\left(\boldsymbol{\zeta}, \Lambda_0\right)
(\boldsymbol{h})=\boldsymbol{\zeta}'\boldsymbol{h_{1}}+\int_{0}^{\tau}h_{2}(u)d\Lambda_0(u)$. 
Then, the parameter $\boldsymbol{\theta}$ can be considered as a functional on $H$, 
and the parameter space $\Theta$ can be considered as a subset of $l^{\infty}(H)$. We 
further define the supremum norm as $\norm{U}_{\infty}=\sup_{\boldsymbol{h} \in H} |
U(\boldsymbol{h})|$ 
on $l^{\infty}(H)$. The score functions $S_{n}$ and $S$ are both maps from $\Theta$ to 
$l^{\infty}(H)$.
 
Next, we verify the four conditions of the Theorem 3.3.1 in \cite{van1996weak}. For verifying the condition (i), let us denote $\boldsymbol{\theta_0}$ 
as $\left(\boldsymbol{\zeta_0},\Lambda_{00}\right)$. We will show that the class
\[
	\mathcal{F}=\left\{l_{1}\left(\boldsymbol{\theta}\right)'\boldsymbol{h_{1}} 
	+ l_{2}\left(\boldsymbol{\theta}\right)[h_{2}] 
	- l_{1}\left(\boldsymbol{\theta_0}\right)'\boldsymbol{h_{1}} 
	- l_{2}\left(\boldsymbol{\theta_0}\right)[h_{2}] 
	| \norm{\boldsymbol{\zeta}-\boldsymbol{\zeta_{0}}}
	+\norm{\Lambda_0-\Lambda_{00}}_{\infty} < \delta, \; \boldsymbol{h} \in H\right\}
\]
is P-Donsker where $\delta>0$ is a constant, by using the Theorem 2.5.6 and 
Theorem 2.7.5 in \cite{van1996weak}.

By the mean value theorem, for any $\boldsymbol{\theta}=(\boldsymbol{\zeta},
\Lambda_0)$, $\boldsymbol{\theta_0}=(\boldsymbol{\zeta_{0}},\Lambda_{00})$ in 
$\Theta=Z \times M$, and $\boldsymbol{h} \in H$, there exists $\boldsymbol{\theta^{*}}=
(1-c)\boldsymbol{\theta} + c \boldsymbol{\theta_0}$ for some $c \in 
\left[0,1\right]$, such that
\begin{equation}\label{eq:15}
\begin{split}
	&l_{1}\left(\boldsymbol{\theta}\right)'\boldsymbol{h_{1}} 
	+ l_{2}\left(\boldsymbol{\theta}\right)[h_{2}] 
	- l_{1}\left(\boldsymbol{\theta_0}\right)'\boldsymbol{h_{1}} 
	- l_{2}\left(\boldsymbol{\theta_0}\right)[h_{2}]\\
	&=(\boldsymbol{\zeta}-\boldsymbol{\zeta_{0}})' 
	\left\{\nabla_{\boldsymbol{\zeta}} 
	m_{1}\left(O | \boldsymbol{\theta^{*}}\right)'\boldsymbol{h_{1}}\right\}
	-(\boldsymbol{\zeta}-\boldsymbol{\zeta_{0}})' 
	\int_{0}^{O}\nabla_{\boldsymbol{\zeta}} 
	m_{2}\left(u | \boldsymbol{\theta^{*}}\right)'\boldsymbol{h_{1}} 
	d\Lambda_0(u)\\
	&-\int_{0}^{O}
	\left\{m_{2}\left(u | \boldsymbol{\theta^{*}}\right)'\boldsymbol{h_{1}} \right\}
	d\left(\Lambda_0-\Lambda_{00}\right)(u)
	-(\boldsymbol{\zeta}-\boldsymbol{\zeta_{0}})' \int_{0}^{O}
	\nabla_{\boldsymbol{\zeta}} m_{3}\left(u | \boldsymbol{\theta^{*}}\right)h_{2}(u)
	d\Lambda_0(u)\\
	&-\int_{0}^{O} 
	\left\{ m_{3}\left(u | \boldsymbol{\theta^{*}}\right)h_{2}(u)\right\}
	d\left(\Lambda_0-\Lambda_{00}\right)(u),
\end{split}
\end{equation}
where $\nabla_{\boldsymbol{\zeta}} m_{1}\left(O | \boldsymbol{\theta^{*}}\right)'
\boldsymbol{h_{1}}$, $\nabla_{\boldsymbol{\zeta}} m_{2}\left(u | 
\boldsymbol{\theta^{*}}\right)'\boldsymbol{h_{1}}$, and $\nabla_{\zeta} m_{3}\left(O | 
\boldsymbol{\theta^{*}}\right)$ denote the derivatives of $m_{1}\left(O | 
\boldsymbol{\theta^{*}}\right)'\boldsymbol{h_{1}}$, $m_{2}\left(u | 
\boldsymbol{\theta^{*}}\right)'\boldsymbol{h_{1}}$, and $m_{3}\left(u | 
\boldsymbol{\theta^{*}}\right)$ with respect to $\boldsymbol{\zeta}$, respectively. 
By the assumptions (C.1) - (C.5), we can show that $|m_{i}\left(u | 
\boldsymbol{\theta^{*}}\right)| \leq c_{i0}e^{c_{i1}\norm{Y}}$ and 
$|\nabla_{\boldsymbol{\zeta}} m_{j}\left(u | \boldsymbol{\theta^{*}}\right)| 
\leq c_{j0}e^{c_{j1}\norm{Y}}$, for constants $c_{ij}$ where $i=1,2,3$ and $j=0,1$. 
Then, from the equation \eqref{eq:15} , we have
\begin{equation}\label{eq:16}
\begin{split}
	&|l_{1}\left(\boldsymbol{\theta}\right)'\boldsymbol{h_{1}} 
	+ l_{2}\left(\boldsymbol{\theta}\right)[h_{2}] 
	- l_{1}\left(\boldsymbol{\theta_0}\right)'\boldsymbol{h_{1}} 
	- l_{2}\left(\boldsymbol{\theta_0}\right)[h_{2}]|\\
	&\leq c_{0}e^{c_{1}\norm{Y}} \left\{\norm{\boldsymbol{\zeta}-\boldsymbol{\zeta_{0}}}
	+|\Lambda_0(O)-\Lambda_{00}(O)|
	+\int_{0}^{\tau}|\Lambda_0(u)-\Lambda_{00}(u)|du\right\}\\
	&+\bigg| \int_{0}^{O} 
	\left\{ m_{3}\left(u | \boldsymbol{\theta^{*}}\right)h_{2}(u)\right\}
	d\left(\Lambda_0-\Lambda_{00}\right)(u)\bigg| 
\end{split}
\end{equation}
for some constants $c_{0}$ and $c_{1}$. For every probability measure $P$, and 
constants $K_{1}$ and $K_{2}$, the entropy number for the class 
$\mathcal{A}=\{\Lambda(\tau)\leq K_{1}\}$ satisfies 
$\log N_{\left[ \cdot \right]}\left(\epsilon,\mathcal{A},L_{2}(P)\right) 
\leq K_{2}/\epsilon$ by the Theorem 2.7.5 in \cite{van1996weak}. 
Hence, to cover the class $\mathcal{A}$, we can find 
$\exp\left\{K_{2}/ \epsilon\right\}$ brackets 
$\left\{\left[L_{i},U_{i}\right]\right\}$ such that $\norm{U_{i}-L_{i}}_{L_{2}(P)} 
\leq \epsilon$ for each $i$. Let $\left\{Q_{1},Q_{2}, \cdots \right\}$ be the 
partition set of $Z$ such that the number of partitions is of the order 
$(1/\epsilon)^{r}$ and for any $\boldsymbol{\zeta}$ and $\boldsymbol{\zeta_{0}}$ 
in the same partition $S_{k}$, their Euclidean distance is less than $\epsilon$. 
Then, $Z \times M$ can be bracket covered by $\left\{Q_{1},Q_{2}, \cdots \right\} 
\times \left\{\left[L_{i},U_{i}\right]\right\}$ and the total number of partition 
is of order $(1/\epsilon)^{r}\exp\left\{1/ \epsilon\right\}$. Hence, from the 
equation \eqref{eq:16}, we obtain $N_{\left[ \cdot \right]}\left(\epsilon,\mathcal{F},
\norm{\cdot}_{L_{2}(P)}\right) \leq O(1)\left(\frac{1}{\epsilon}\right)^{r}
\exp\left\{1/\epsilon\right\}$. Hence, $\mathcal{F}$ is P-Donsker by the 
Theorem 2.5.6 in \cite{van1996weak}.

Moreover, it follows that the right-hand side 
of the equation \eqref{eq:16} uniformly converges to $0$ as $\norm{\boldsymbol{\zeta}-
\boldsymbol{\zeta_{0}}}+\norm{\Lambda_0-\Lambda_{00}}_{\infty} \rightarrow 0$. 
That is, 
\[
	\sup_{\boldsymbol{h} \in H} E_{0}\left\{
	l_{1}\left(\boldsymbol{\theta}\right)'\boldsymbol{h_{1}} 
	+ l_{2}\left(\boldsymbol{\theta}\right)[h_{2}] 
	- l_{1}\left(\boldsymbol{\theta_0}\right)'\boldsymbol{h_{1}} 
	- l_{2}\left(\boldsymbol{\theta_0}\right)[h_{2}]
	\right\}^2 \rightarrow 0 
\]
as $\norm{\boldsymbol{\zeta}-\boldsymbol{\zeta_{0}}}+\norm{\Lambda_0-
\Lambda_{00}}_{\infty} \rightarrow 0$. Thus, the condition (i) is satisfied by the
Lemma 3.3.5 of \cite{van1996weak}.

The condition (ii) in the Theorem 3.3.1 of \cite{van1996weak}
follows from the fact that $\sqrt{n}\left(S_{n}-S\right)(\boldsymbol{\theta_{0}})$ 
converges in distribution to a tight Gaussian process $\mathcal{G}$ in $l^{\infty}(H)$, 
which is a direct conclusion by combining the Donsker theorem and the verification of 
condition (i).

To verify the condition (iii), we first notice that the derivative of the score 
function $S\left(\boldsymbol{\theta_{0}}\right)(\boldsymbol{h})$ is given by
\[
   D_{\boldsymbol{\theta}}S\left(\boldsymbol{\theta_{0}}\right)(\boldsymbol{h})
   =\frac{\partial}{\partial t}S\left(\boldsymbol{\theta_{0}}+
   t\boldsymbol{\theta}\right)(\boldsymbol{h})\bigg|_{t=0},
\]
where $\boldsymbol{\theta_{0}}+t\boldsymbol{\theta}=\left(\boldsymbol{\zeta_{0}}+
t\boldsymbol{\zeta},\Lambda_{00}(u)+t\Lambda_0(u)\right)$.
Then, it can be calculated that
\begin{equation}
\begin{split}
	&\frac{\partial}{\partial t}S\left(\boldsymbol{\theta_{0}}+
        t\boldsymbol{\theta}\right)(\boldsymbol{h})\bigg|_{t=0}\\
	&=E_{0}\left\{\nabla_{\boldsymbol{\zeta}} 
	m_{1}\left(O | \boldsymbol{\theta_{0}}\right)'\boldsymbol{h_{1}}
	-\int_{0}^{O}\nabla_{\boldsymbol{\zeta}} 
	m_{2}\left(u | \boldsymbol{\theta_{0}}\right)'\boldsymbol{h_{1}} 
	d\Lambda_{00}(u)\right.\\
	&\left. -\int_{0}^{\tau}I(u \leq O)\nabla_{\boldsymbol{\zeta}} 
	m_{3}\left(u | \boldsymbol{\theta_{0}}\right)h_{2}(u)
	d\Lambda_{00}(u)\right\}\\
	&-\int_{0}^{\tau}E_{0}\left[I(u \leq O)\left\{
	m_{2}\left(u | \boldsymbol{\theta_{0}}\right)'\boldsymbol{h_{1}}
	+m_{3}\left(u | \boldsymbol{\theta_{0}}\right)h_{2}(u)\right\}\right]
	d\Lambda_{00}(u),\nonumber
\end{split}
\end{equation}
where $\nabla_{\boldsymbol{\zeta}} m_{1}\left(O | \boldsymbol{\theta_{0}}\right)'
\boldsymbol{h_{1}}$, $\nabla_{\boldsymbol{\zeta}} m_{2}\left(u | 
\boldsymbol{\theta_{0}}\right)'\boldsymbol{h_{1}}$, and $\nabla_{\zeta} m_{3}\left(O 
| \boldsymbol{\theta_{0}}\right)$ denote the derivatives of $m_{1}\left(O | 
\boldsymbol{\theta_{0}}\right)'\boldsymbol{h_{1}}$, $m_{2}\left(u | 
\boldsymbol{\theta_{0}}\right)'\boldsymbol{h_{1}}$, and $m_{3}\left(u | 
\boldsymbol{\theta_{0}}\right)$ with respect to $\boldsymbol{\zeta}$, respectively. 
By the first-order Taylor's expansion of the score function 
$S\left(\boldsymbol{\theta_{0}}+\epsilon\boldsymbol{\theta}\right)(\boldsymbol{h})$ 
around $\theta_{0}$, we get
\[
	S\left(\boldsymbol{\theta_{0}}+\epsilon\boldsymbol{\theta}\right)(\boldsymbol{h})
	=S\left(\boldsymbol{\theta_{0}}\right)(\boldsymbol{h}) 
	+ D_{\epsilon\boldsymbol{\theta}}
	S\left(\boldsymbol{\theta_{0}}\right)(\boldsymbol{h})
	+\littleo\left(\epsilon\right).
\]
Then, 
\[
	\frac{\norm{S\left(\boldsymbol{\theta_{0}}+
        \epsilon\boldsymbol{\theta}\right)(\boldsymbol{h})
	-S\left(\boldsymbol{\theta_{0}}\right)(\boldsymbol{h})
	-D_{\epsilon\boldsymbol{\theta}}
	S\left(\boldsymbol{\theta_{0}}\right)(\boldsymbol{h})}_{\infty}}
	{\epsilon}
\]
converges uniformly to $0$ with respect to $\boldsymbol{\theta}$ when it takes
any element of bounded subsets of all finite linear combinations of elements of 
$\Theta$, as $\epsilon \rightarrow 0$. Hence, $S(\boldsymbol{\theta})$ is 
Fr\'echet-differentiable at $\boldsymbol{\theta_{0}}$ and the Fr\'echet 
derivative is $\nabla{S}_{\boldsymbol{\theta_{0}}}=D_{\boldsymbol{\theta}}
S\left(\boldsymbol{\theta_{0}}\right)$. Then, $\nabla{S}_{\boldsymbol{\theta_{0}}}
\left(\boldsymbol{\theta}-\boldsymbol{\theta_{0}}\right)(\boldsymbol{h})$ can be 
written as
\[
   (\boldsymbol{\zeta}-\boldsymbol{\zeta_{0}})'\omega_{1}(\boldsymbol{h_{1}},h_{2})
   +\int_{0}^{\tau}\omega_{2}(\boldsymbol{h_{1}},h_{2})d(\Lambda_0-\Lambda_{00})(u),
\]
where
\begin{equation}\label{eq:17}
\begin{split}
	&\omega_{1}(\boldsymbol{h_{1}},h_{2})=
	E_{0}\left\{\nabla_{\boldsymbol{\zeta}} 
	m_{1}\left(O | \boldsymbol{\theta_0}\right)'\boldsymbol{h_{1}}
	-\int_{0}^{O}\nabla_{\boldsymbol{\zeta}} 
	m_{2}\left(u | \boldsymbol{\theta_0}\right)'\boldsymbol{h_{1}} d\Lambda_{00}(u)\right.\\
	&\qquad \qquad \left.-\int_{0}^{\tau}I(u \leq O)\nabla_{\boldsymbol{\zeta}} 
	m_{3}\left(u | \boldsymbol{\theta_0}\right)h_{2}(u)d\Lambda_{00}(u)\right\}\\
	&\omega_{2}(\boldsymbol{h_{1}},h_{2})
	=-E_{0}\left[I(u \leq O)\left\{m_{2}\left(u | 
        \boldsymbol{\theta_0}\right)'\boldsymbol{h_{1}}
	+m_{3}\left(u | \boldsymbol{\theta_0}\right)h_{2}(u)\right\}\right].
\end{split}
\end{equation}
Note that the linear operator $\omega(\boldsymbol{h})=\left(\omega_{1}
(\boldsymbol{h_{1}},h_{2}),\omega_{2}(\boldsymbol{h_{1}},h_{2})\right)$ is a bounded 
operator from $H \subset R^{r} \times \boldsymbol{BV}_{\left[0,\tau\right]}$ to 
$R^{r} \times \boldsymbol{BV}_{\left[0,\tau\right]}$ where 
$\boldsymbol{BV}_{\left[0,\tau\right]}$ contains all the functions with finite 
total variation in $\left[0,\tau\right]$. Note also that 
$\nabla{S}_{\boldsymbol{\theta_{0}}}$ is continuous invertible if 
\[
   \inf_{\boldsymbol{\theta} \in \Theta} \frac{\norm{\nabla{S}_{\boldsymbol{\theta_{0}}}
   \left(\boldsymbol{\theta}-\boldsymbol{\theta_{0}}\right)(\boldsymbol{h})}_{\infty}}
   {\norm{\left(\boldsymbol{\theta}-\boldsymbol{\theta_{0}}\right)
   (\boldsymbol{h})}_{\infty}} \geq \epsilon, \qquad 
   \mbox{for some } \epsilon > 0.
\]
Since $\left(\boldsymbol{\theta}-\boldsymbol{\theta_{0}}\right)(\boldsymbol{h})$ 
belongs to $l^{\infty}(H)$, we can regard $\nabla{S}_{\boldsymbol{\theta_{0}}}
\left(\boldsymbol{\theta}-\boldsymbol{\theta_{0}}\right)(\boldsymbol{h})$ as an 
operator from $l^{\infty}(H)$ to $l^{\infty}(H)$.
If $\omega(\boldsymbol{h})=\left(\omega_{1}(\boldsymbol{h_{1}},h_{2}),
\omega_{2}(\boldsymbol{h_{1}},h_{2})\right)$ is a continuous invertible operator, 
which means that there exists a constant $\epsilon > 0$ such that 
$\omega^{-1}\left(\epsilon H\right) \subset H$, then for any 
$\left(\boldsymbol{\theta}-\boldsymbol{\theta_{0}} \right)(\boldsymbol{h}) \in l^{\infty}(H)$, we have
\begin{equation}
\begin{split}
	&\inf_{\boldsymbol{\theta} \in \Theta} 
	\frac{\norm{\nabla{S}_{\boldsymbol{\theta_0}}\left(\boldsymbol{\theta}-
        \boldsymbol{\theta_0}\right)(\boldsymbol{h})}_{\infty}}
	{\norm{\left(\boldsymbol{\theta}-\boldsymbol{\theta_0}\right)
        (\boldsymbol{h})}_{\infty}}\\
	&=\inf_{\boldsymbol{\theta} \in \Theta} 
	\frac{\sup_{\boldsymbol{h} \in H}{|\nabla{S}_{\boldsymbol{\theta_0}}
        \left(\boldsymbol{\theta}-\boldsymbol{\theta_0}\right)(\boldsymbol{h})|}}
	{\norm{\left(\boldsymbol{\theta}-\boldsymbol{\theta_0}\right)
        (\boldsymbol{h})}_{\infty}}\\
	&=\inf_{\boldsymbol{\theta} \in \Theta} 
	\frac{\sup_{\boldsymbol{h} \in H}{|(\boldsymbol{\zeta}-\boldsymbol{\zeta_{0}})'
        \omega_{1}(\boldsymbol{h_{1}},h_{2})+\int_{0}^{\tau}
        \omega_{2}(\boldsymbol{h_{1}},h_{2})d(\Lambda_0-\Lambda_{00})(u)|}}
	{\norm{\left(\boldsymbol{\theta}-\boldsymbol{\theta_0}\right)
        (\boldsymbol{h})}_{\infty}}\\
	&\geq \inf_{\boldsymbol{\theta}\in \Theta} 
	\frac{\epsilon{\norm{\left(\boldsymbol{\theta}-\boldsymbol{\theta_0}\right)
        (\boldsymbol{h})}}_{\infty}}
	{\norm{\left(\boldsymbol{\theta}-\boldsymbol{\theta_0}\right)
        (\boldsymbol{h})}_{\infty}}
	= \epsilon.\nonumber
\end{split}
\end{equation}
Thus, showing that $\nabla{S}_{\boldsymbol{\theta_{0}}}$ is a continuous invertible
operator is equivalent to showing that $\omega(\boldsymbol{h})$ is a continuous 
invertible operator. By the results on page 424 of \cite{van1998asymptotic}, showing that
$\omega(\boldsymbol{h})$ is a continuous invertible operator is equivalent to showing
that $\omega(\boldsymbol{h})$ is one-to-one on its range and it can be written as 
the sum of a continuously invertible operator and a compact linear operator.
To prove that $\omega(\boldsymbol{h})$ is one-to-one, we need to show that if 
$\omega(\boldsymbol{h})=0$ for some $\boldsymbol{h} \in H$, then $\boldsymbol{h}=0$. 
Assume that $\omega(\boldsymbol{h_1},h_{2})=0$ for some $\boldsymbol{h}=
(\boldsymbol{h_1},h_{2}) \in H$. Then, 
\begin{equation}\label{eq:18}
	\boldsymbol{h_1}'\omega_{1}(\boldsymbol{h_1},h_{2})
	+\int_{0}^{\tau}h_{2}(u)\omega_{2}(\boldsymbol{h_{1}},h_{2})d\Lambda_{0}(u)=0
\end{equation}
and both $\omega_{1}(\boldsymbol{h_{1}},h_{2})$ and $\omega_{2}(\boldsymbol{h_{1}},h_{2})$ 
are identically $0$. The equation \eqref{eq:18} can be re-expressed as 
$\nabla{S}_{\theta_{0}}\left(\boldsymbol{h_{1}},\int_{0}^{\tau}h_{2}(u)d\Lambda_{0}(u)\right)
(\boldsymbol{h})=0$, and the left-hand-side is the negative information matrix. 
Thus, it follows that the score function $l_{1}\left(\boldsymbol{\theta_0}\right)'
\boldsymbol{h_{1}} + l_{2}\left(\boldsymbol{\theta_0}\right)[h_{2}]$ is almost everywhere 
$0$. By using the same arguments as those in the proof of Lemma 2, we can obtain that 
$\boldsymbol{h_{1}}=0$ and $h_{2}=0$.

Next, we need to show that the operator $\omega$ can be written as the sum of a 
continuously invertible operator $A$ and a compact linear operator 
$B=\left(B_{1},B_{2}\right)$; that is, $\omega\left(\boldsymbol{h_{1}},h_{2}\right)
=A\left(\boldsymbol{h_{1}},h_{2}\right)+\left(B_{1}\left(\boldsymbol{h_{1}},h_{2}\right),
B_{2}\left(\boldsymbol{h_{1}},h_{2}\right)\right)$, where
\begin{equation}
\begin{split}
	A\left(\boldsymbol{h_{1}},h_{2}\right)
	&=\left(\boldsymbol{h_{1}}, -P_{0}\left[I(u \leq O)m_{3}\left(u | 
        \boldsymbol{\theta_0}\right)h_{2}(u)\right]\right) \\
	B_{1}\left(\boldsymbol{h_{1}},h_{2}\right)
	&=\omega_{1}(\boldsymbol{h_{1}},h_{2})-\boldsymbol{h_{1}}\\
	B_{2}\left(\boldsymbol{h_{1}},h_{2}\right)
	&=\omega_{2}(\boldsymbol{h_{1}},h_{2}) + P_{0}\left[I(u \leq O)m_{3}\left(u | 
        \boldsymbol{\theta_0}\right)h_{2}(u)\right]. \nonumber
\end{split}
\end{equation}
Note that the linear operator $\omega(\boldsymbol{h})=\left(\omega_{1}
(\boldsymbol{h_{1}},h_{2}),\omega_{2}(\boldsymbol{h_{1}},h_{2})\right)$ is a bounded 
operator from $H \subset R^{r} \times \boldsymbol{BV}_{\left[0,\tau\right]}$ to $R^{r} \times 
\boldsymbol{BV}_{\left[0,\tau\right]}$. By the assumptions (C.1) -(C.5), it follows that 
$A^{-1}$ is a bounded linear operator. Thus, $A$ is continuously invertible. Because 
a bounded linear operator with finite-dimensional range is compact, $B_{1}$ is compact. 
By the Arzel\'a-Ascoli theorem, $B_{2}$ is a compact operator from $R^{r} \times 
\boldsymbol{BV}_{\left[0,\tau\right]}$ to $\boldsymbol{BV}_{\left[0,\tau\right]}$ since 
the range of $B_{2}$ is a continuously 
differentiable function in $\left[0,\tau\right]$. Hence, it follows that $B$ is 
a compact operator. Thus, we 
have proved that operator $\omega$ is the sum of a continuously invertible operator 
and a compact linear operator, and we conclude that  $\omega(\boldsymbol{h})$ is 
a continuous invertible operator. It follows directly that 
$\nabla{S}_{\boldsymbol{\theta_0}}$ is continuous invertible.

The condition (iv) in the Theorem 3.3.1 of \cite{van1996weak} can
be verified by the following two results. First, from Lemma 1 about
the existence of MLEs and the preceding arguments in this proof, we 
can conclude that $S_{n}(\hat{\boldsymbol{\theta}})\left(\boldsymbol{h}\right)=0$. 
Second, from the proof of Lemma 2, we can easily obtain the result that 
$S\left(\boldsymbol{\theta_0}\right)\left(\boldsymbol{h}\right)=0$. Thus, the
condition (iv) is satisfied.

So far, we have verified all four conditions in the Theorem 3.3.1 of \cite{van1996weak}. From that theorem, we have the conclusion that
\begin{equation}\label{eq:19}
	- \nabla{S}_{\boldsymbol{\theta_{0}}}
	\sqrt{n}(\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta}, 
	\hat{\Lambda}_0-\Lambda_0)(\boldsymbol{h})
	=\sqrt{n}\left(S_{n}(\boldsymbol{\theta_{0}})-S(\boldsymbol{\theta_{0}})\right)(\boldsymbol{h}) + o_{P}(1)
\end{equation}
where $o_{P}(1)$ denotes convergence to zero in probability in the space $l^{\infty}(H)$,
and
\begin{equation}\label{eq:20}
	- \nabla{S}_{\boldsymbol{\theta_{0}}}
	\sqrt{n}(\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta}, 
	\hat{\Lambda}_0-\Lambda_0)(\boldsymbol{h})
	=\sqrt{n}(\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta})'\omega_{1}(\boldsymbol{h}) 
	+\int_{0}^{\tau}\omega_{2}(\boldsymbol{h})\sqrt{n}d(\hat{\Lambda}_0-\Lambda_0)(u).
\end{equation}
Hence, we can conclude that 
$\sqrt{n}(\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta}, 
\hat{\Lambda}_0-\Lambda_0)(\boldsymbol{h})$ converges in probability to 
$- \nabla{S}_{\boldsymbol{\theta_{0}}}^{-1}\mathcal{G}(\boldsymbol{h})$ where $\mathcal{G}(\boldsymbol{h})$ is the limiting distribution of $\sqrt{n}\left(S_{n}(\boldsymbol{\theta_{0}})-S(\boldsymbol{\theta_{0}})\right)(\boldsymbol{h})$. Moreover, since $\omega(\boldsymbol{h})$ is a continuous 
invertible operator, there exists $\boldsymbol{g}=(\boldsymbol{g}_1,g_2)$ in $H$ 
such that $\boldsymbol{h}=\omega^{-1}(\boldsymbol{g})$. Then, from the equation \eqref{eq:19}, 
we have
\begin{equation}\label{eq:21}
	- \nabla{S}_{\boldsymbol{\theta_{0}}}
	\sqrt{n}(\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta}, 
	\hat{\Lambda}_0-\Lambda_0)(\boldsymbol{h})
	=\sqrt{n}\left(S_{n}(\boldsymbol{\theta_{0}})-S(\boldsymbol{\theta_{0}})\right)
        \left(\omega^{-1}(\boldsymbol{g})\right) + o_{P}(1).
\end{equation}
The right hand side of the equation \eqref{eq:21} converges to 
$\mathcal{G}(\omega^{-1}(\boldsymbol{g}))$, which has mean zero and variance 
$\omega_{1}^{-1}(\boldsymbol{g})'\boldsymbol{g}_{1} +
\int_{0}^{\tau}\omega_{2}^{-1}(\boldsymbol{g})g_{2}(u)d\Lambda_0(u)$. Moreover, 
from the equation \eqref{eq:20}, we have
\begin{equation}\label{eq:22}
   - \nabla{S}_{\boldsymbol{\theta_{0}}}\sqrt{n}(\hat{\boldsymbol{\zeta}}-
   \boldsymbol{\zeta}, \hat{\Lambda}_0-\Lambda_0)(\boldsymbol{h})=
   \sqrt{n}(\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta})'\boldsymbol{g}_{1} 
   +\int_{0}^{\tau}g_{2}(u)\sqrt{n}d(\hat{\Lambda}_0-\Lambda_0)(u),
\end{equation}
which implies that $\sqrt{n}(\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta},
\hat{\Lambda}_0-\Lambda_0)(\boldsymbol{g})$ converges in probability to 
$- \nabla{S}_{\boldsymbol{\theta_{0}}}^{-1}\mathcal{G}(\boldsymbol{g})$. Thus, by the 
equations \eqref{eq:21} and \eqref{eq:22}, we have 
$- \nabla{S}_{\boldsymbol{\theta_{0}}}^{-1}\mathcal{G}(\boldsymbol{g})=
\mathcal{G}(\omega^{-1}(\boldsymbol{g}))$. Hence, we can conclude that 
$\sqrt{n}(\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta}, 
\hat{\Lambda}_0-\Lambda_0)$ converges weakly to a Gaussian process G in 
$l^{\infty}\left(H\right)$ with mean zero and the covariance process
\[
   cov\left[G(\boldsymbol{g}),G(\boldsymbol{g}^{*})\right]=
   \omega_{1}^{-1}(\boldsymbol{g}^{*})'\boldsymbol{g}_{1} +
   \int_{0}^{\tau}\omega_{2}^{-1}(\boldsymbol{g}^{*})g_{2}(u)d\Lambda_0(u),
\]
where $(\omega_{1}^{-1},\omega_{2}^{-1})$ is the inverse of the continuously 
invertible linear operator $\left(\omega_{1},\omega_{2}\right)$. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\section*{References}
\bibliographystyle{tfnlm}
\bibliography{ref2}

\end{document}
