Nonparametric inference for continuous-time event counting and link-based dynamic network models

TL;DR

This paper introduces a nonparametric likelihood-based method for modeling and estimating dynamic event counting and link-based networks using kernel smoothing, with theoretical analysis and empirical validation.

Contribution

It develops a nonparametric estimation framework for continuous-time network models and analyzes their asymptotic properties as the network size grows.

Findings

Asymptotic properties of estimators are rigorously established.

Finite sample performance demonstrated on bike share data.

Method effectively captures dynamic network behaviors.

Abstract

A flexible approach for modeling both dynamic event counting and dynamic link-based networks based on counting processes is proposed, and estimation in these models is studied. We consider nonparametric likelihood based estimation of parameter functions via kernel smoothing. The asymptotic behavior of these estimators is rigorously analyzed by allowing the number of nodes to tend to infinity. The finite sample performance of the estimators is illustrated through an empirical analysis of bike share data.

Figures

Equations

$$N_{n,ij}(t)=\#\{\text{interaction events between }i\text{ and }j\text{ up to and including time }t\}.$$

$${\lambda }_{n,ij}(\theta ,t):=\Phi (\theta (t);{\left(X_{n,ij}(s)\right)}_{i,j=1,...,n}:s\le t),$$

$$\Phi (\theta ;{\left(C_{n,ij}(s),X_{n,ij}(s)\right)}_{i,j=1,...,n}:s\le t)=C_{n,ij}(t)\exp ({\theta }^T(t)X_{n,ij}(t)).$$

$${\ell }_T(\theta )=\sum _{0<t\le T}\sum _{(i,j)\in L_n}\Delta N_{n,ij}(t)\log {\lambda }_{n,ij}(\theta ,t)-{\int }_0^T\sum _{(i,j)\in L_n}{\lambda }_{n,ij}(\theta ,t)dt,$$

$$\hat{\theta }:=\arg \underset{\theta \in \Theta }{\max }{\ell }_T(\theta ),$$

$${\ell }_T(\mu ,t_0)$$

$$-{\int }_0^T\sum _{(i,j)\in L_n}\frac{1}{h}K\left(\frac{t-t_0}{h}\right){\lambda }_{n,ij}(\mu ,t)dt,$$

$$\hat{\theta }(t_0)=\arg \underset{\theta \in \Theta }{\max }{\ell }_T(\theta ,t_0),$$

$$\displaystyle\lambda_{n,ij}(\theta,t)=C_{n,ij}(t)\exp\big{\{}\theta(t)^{T}X_{n,ij}(t)\big{\}}.$$

$${\ell }_T(\theta ,t_0)=\sum _{i,j=1}^n{\int }_0^T\frac{1}{h}K\left(\frac{t-t_0}{h}\right){\theta }^T{X}_{n,ij}(t)d{N}_{n,ij}(t)$$

$$-\sum _{i,j=1}^n{\int }_0^T\frac{1}{h}K\left(\frac{t-t_0}{h}\right)C_{n,ij}(t)\exp ({\theta }^T{X}_{n,ij}(t))dt.$$

$$\int \underset{s\in U}{\sup }\{\left(1+\Vert y\Vert +\Vert y{\Vert }^2+\Vert y{\Vert }^3\right)\left|f_s(y)\right|+\left(1+\Vert y\Vert +\Vert y{\Vert }^2\right)\left|{\partial }_s{f}_s(y)\right|$$

$$+\left(1+\Vert y\Vert \right)\left|\partial s^2{f}_s(y)\right|+\Vert y{\Vert }^2\cdot f_{s,t_0}(y)\}\cdot \exp (\tau \cdot \Vert y\Vert )\mathrm{d}\mu (y)<\infty ,$$

$$\displaystyle\sup_{s,t\in U_{h}}\mathbb{E}\bigg{(}\|X_{n,ij}(s)\|^{2}\cdot\|X_{n,kl}(t)\|^{2}$$

$$\cdot e^{\tau (\Vert X_{n,ij}(s)\Vert +\Vert X_{n,kl}(t)\Vert )}|C_{n,ij}(t_0)=1,C_{n,kl}(t_0)=1)=O(1).$$

$$\displaystyle\sup_{s\in U_{h}}\mathbb{E}\left(\|X_{n,12}(s)\|^{k}e^{\tau\|X_{n,12}(s)\|}\bigg{|}C_{n,12}(s)=1,C_{n,12}(t_{0})=0\right)=O(1),$$

$$\displaystyle\mathbb{E}\bigg{(}\sup_{s\in U_{h}}\left[\|X_{n,12}(s)\|+\|X_{n,12}(s)\|^{2}+\|X_{n,12}(s)\|^{3}+\|X_{n,12}(s)\|^{4}\right]$$

$$\displaystyle\hskip 142.26378pt\cdot e^{\tau\|X_{n,12}(s)\|}\bigg{|}C_{n,12}(s)=1\bigg{)}<+\infty.$$

$${\int }_{-1}^{1}w(u)\frac{\mathbb{P}(C_{n,12}(t_0+uh)=1)}{\mathbb{P}(C_{n,12}(t_0)=1)}\mathrm{d}u\to 1\text{ as }n\to \infty .$$

$$A_{n,ij,kl}$$

$$A_{n,ij,kl}=\left\{\begin{array}[]{ccc}o(n^{2})&\text{for}&|\{i,j\}\cap\{k,l\}|=2,\\ o(n)&\text{for}&|\{i,j\}\cap\{k,l\}|=1,\\ 1+o(1)&\text{for}&|\{i,j\}\cap\{k,l\}|=0.\end{array}\right.$$

$${\int }_0^T\frac{1}{h}K\left(\frac{s-t_0}{h}\right)\frac{\mathbb{P}(C_{n,12}(t_0)=0,C_{n,12}(s)=1)}{\mathbb{P}(C_{n,12}(t_0)=1)}ds=O(h),$$

$$\frac{\mathbb{P}(C_{n,ij}(t_{0})=1,C_{n,kl}(t_{0})=1)}{\big{(}\mathbb{P}(C_{n,12}(t_{0})=1)\big{)}^{2}}=O(1).\\$$

$$\frac{\mathbb{P}(C_{n,ij}(t_0)=0,C_{n,kl}(s)=1)}{\mathbb{P}(C_{n,12}(t_0)=1)}=\frac{\mathbb{P}({\delta }_{n,ij}\in [s-{\delta }_0,t_0-{\delta }_0))}{\mathbb{P}({\delta }_{n,ij}\ge t_0-{\delta }_0)}.$$

$$\mathbb{P}(C_{n,12}=1)$$

$$\mathbb{P}(C_{n,12}=1,C_{n,23}=1)$$

$$\mathbb{P}(C_{n,12}=1,C_{n,34}=1)$$

$$\theta \mapsto {\int }_0^T\frac{1}{h}K\left(\frac{s-t_0}{h}\right)g(\theta ,s)ds,$$

$${\tau }_{n,ij}(\theta ,s)$$

$$g(\theta ,t)$$

$$=$$

$$f_{n,1}(\theta ,s,t|(i,j),(k,l))$$

$$f_2(\theta ,t)$$

$$r_{n,ij}^{(a)}(s)$$

$$f_1({\theta }_0(t_0),0)=f_2({\theta }_0(t_0),t_0{)}^2.$$

$$\iint\limits_{[-1,1]^{2}}K(u)K(v)\mathbb{E}\left(\rho_{n,ijkl}^{(a)}(u,v)|C_{n,ij}(t_{0})=1,C_{n,kl}(t_{0})=1\right)\mathrm{d}u\mathrm{d}v=o\big{(}\big{(}l_{n}h\big{)}^{-1}\big{)}.$$

$$\mathbb{E}(A_{i,j}{A}_{k,l})$$

Full text

Nonparametric inference for continuous-time event counting and link-based dynamic network models

Alexander Kreiß

University of Mannheim

Department of Economics

L7, 3-5

68161 Mannheim

Germany

kreiss@uni-mannheim.de

  

Enno Mammen

Heidelberg University

Institute for Appplied Mathematics

Im Neuenheimer Feld 205

69120 Heidelberg

Germany

mammen@math.uni-heidelberg.de

  

Wolfgang Polonik

Department of Statistics

University of California, Davis

One Shields Ave.

Davis, CA 95616

USA

wpolonik@ucdavis.edu

A flexible approach for modeling both dynamic event counting and dynamic link-based networks based on counting processes is proposed, and estimation in these models is studied. We consider nonparametric likelihood based estimation of parameter functions via kernel smoothing. The asymptotic behavior of these estimators is rigorously analyzed in an asymptotic framework where the number of nodes tends to infinity. The finite sample performance of the estimators is illustrated through an empirical analysis of bike share data.

1 Introduction

In this paper we present a modeling approach that can be applied to both dynamic interactions in networks as well as dynamic link deletion and addition in networks. The case of dynamic interactions considers a network as a collection of actors who can cause instantaneous interactions. Both directed and undirected interactions are considered. In the model we assume that the time at which an interaction happens, and the pair of actors involved in this interaction, are random. We call this model a dynamic network interaction model. The Enron e-mail data set provides a typical example for a data set that can be modeled in such a way. Here one person sending an e-mail to another person is interpreted as the interaction. While such interactions can be thought of as edges between two nodes, and while the nodes themselves persist over time, each such edge only exists for an infinitesimal time. In contrast to interaction events, we also consider models for connections between the actors of a network that persist over a longer time period. In this case, connections comprise four quantities: sender, receiver (a pair of actors), and the lifetime of each connection (start and end). Here we speak of dynamic networks. Examples for such models are social networks with links indicating an ongoing friendship between two actors. The notions of dynamic network interactions and of dynamic networks are different but they are closely related. A dynamic network defines two network interaction models, one given by the starting time and one given by the ending time of a connection. Furthermore, a network interaction model of e-mails defines a dynamic network by aggregation, where the network shows a connection between two actors as long as they have exchanged e-mails over a certain time period in the past.

In this paper we mainly consider dynamic network interaction models, because dynamic networks can be understood as two dynamic network interaction models. In the model the distribution of the next event conditional on the past is expressed as depending on two quantities: The covariates and the parameter function. The covariates are random but observed. They summarize the relevant history of past interactions. Naturally they are functions of the time. How the covariates influence the distribution is regulated by an unobserved, deterministic parameter function. Thus, the model allows for the effect of the covariates to be changing over time. It is this effect (the parameter function) that is our primary interest of estimation.

In this paper we will make use of counting process models from survival analysis that we adapt for our network models. We will develop asymptotic theory for a kernel-based estimator of parameter functions in the models. The estimator is based on a localized likelihood criterion.

1.1 Literature Review and Related Work

Random networks/graphs have been considered in various scientific branches since a long time, in particular in the social sciences (c.f. the textbook [37]). The importance of the analysis of random networks within the more statistical and machine learning literature is more recent, but corresponding literature is significant by now (e.g. see [8, 18]). The reason for this increase in significance is not least due to the development of modern technologies that lead to the ever increasing number of complex data sets that are encoding relational structures. Examples for real network data can be found at the Koblenz Network Collection KONECT, the European network data collection SocioPatterns, the MIT based collection Reality Commons, the data sets made available by the Max Planck Institute for Softwaresysteme (MPI-SWS), or the Stanford Large Network Dataset Collection (SNAP). In this paper, we will use the Capital Bikeshare Performance Data (see http://www.capitalbikeshare.com/system-data), which is a data set collected on the Washington, DC bikeshare system.

Modeling and analyzing dynamic random networks is challenging as networks can have a multitude of different topological properties. In the literature such topological properties are measured by various quantities, including the flow through the network, the degree distribution, centrality, the existence of hubs, sparsity etc. Time-varying or dynamic random networks appear quite naturally, even though the dynamic aspect significantly adds to the complexity of modeling and analyzing the networks. Early work on networks involving temporal structures can already be found in [16], who consider a discrete time Markov process for friendships (links). Other relevant literature using discrete time settings include work on dynamic exponential random graph models [30, 31, 9, 1, 10, 22, 23], dynamic infinite relational models [14], dynamic block models [12, 38, 39, 40], dynamic nodal states models [20, 19], various dynamic latent features model [7], dynamic multi-group membership models [17], dynamic latent space models [6, 32], and dynamic Gaussian graphical models [42, 21]. Also time-continuous models have been discussed in the literature. They include link-based continuous-time Markov processes [36, 25, 24], actor-based continuous-time Markov processes [33, 34, 35], and also models based on counting processes as in [29], who are considering the modeling of network interaction data, and [4], who applies such a model to radio communication data. Link Prediction, a problem related to the analysis of dynamic networks, received quite some attention in the computer science literature (e.g. see [26, 3]).

1.2 Our Work

In this paper, we study a network model under the asymptotics that the network size $n$ (the number of actors) is growing to infinity. So let $V_{n}:=\{1,...,n\}$ denote the set of $n$ actors (also called agents or nodes), and let $L_{n}\subseteq\{(i,j):i<j,\,i,j\in V_{n}\}$ denote the set of all possible links among them. For directed networks $L_{n}$ is the set of all ordered pairs. Furthermore, let $G_{n}:=(V_{n},L_{n})$ be the corresponding graph. For each pair of actors $(i,j)\in L_{n}$ we denote by $N_{n,ij}\colon[0,T]\to\mathbb{N}$ the number of interactions between these two

[TABLE]

We assume that for $(i,j)\in L_{n}$, the processes $N_{n,ij}$ are one-dimensional counting processes with respect to an increasing, right continuous, complete filtration ${\cal F}_{t}$, $t\in[0,T]$, i.e. the filtration obeys the conditions habituelles, see [2], pp. 60. The $\sigma$-field ${\cal F}_{t}$ contains all information available up to the time point $t$. For simplicity we formulate all results for undirected interactions only, i.e, we assume that $N_{n,ij}=N_{n,ji}$ for all pairs $(i,j)\in L_{n}$. All results can be formulated for the directed case as well (see also the discussion in Section 2.3 after assumption (A2)).

In our approach, we model the intensities of the counting processes $N_{n,ij}$ at time $t$ only for a random subset of edges $\{(i,j)\in L_{n}:C_{n,ij}(t)=1\}$. The functions $C_{n,ij}(t)\in\{0,1\}$ are indicator functions assumed to be predictable with respect to the filtration ${\cal F}_{t}$, and they determine the active part of the network. Our aim is to model the active part. The definition of the active part depends on the application. For instance, modeling the edges between actors $i$ and $j$ might depend on whether $i$ and $j$ have had a low or a high interaction intensity in the past. We will come back to this point later.

For the set $\{(i,j)\in L_{n}:C_{n,ij}(t)=1\}$ the intensities of the counting processes $N_{n,ij}$ are modeled by

[TABLE]

where $X_{n,ij}\colon[0,T]\to\mathbb{R}^{q}$ are $\mathcal{F}_{t}$-predictable covariates, $\Phi$ is a link function, and the parameter function $\theta\colon[0,T]\to\mathbb{R}^{q}$ is the target of our estimation method. The presented above approach is quite flexible and general. In order to be more specific, and for modeling reasons explained in Section 2.1, we will, in the following, assume that $\Phi$ has the following Cox-type form

[TABLE]

Butts suggests in [4] the same modeling framework with a constant parameter in an empirical analysis to radio communication data from responders to the World Trade Center Disaster. In this paper we extend the work of [4] to time varying parameter functions and provide asymptotic theory. Another related model can be found in [25]. Perry and Wolfe studied in [29] a model similar to (1.1) with constant parameters. In their specification the intensity was equal to $\lambda(t)\exp(\theta^{T}X_{n,ij}(t))$ where $\lambda(t)$ is an unknown baseline hazard. They developed asymptotic theory for maximum partial likelihood estimators of $\theta$ in an asymptotic framework where the time horizon $T$ converges to infinity. Our work was motivated by their research but it differs in several respects. First of all we allow the parameters to change over time. Our estimates of the parameter functions can be used for statistical inference on time changes in the effects of covariates. Furthermore, by choosing the first component $X_{n,ij}^{(1)}=1$ our model includes the time-varying baseline intensity $e^{\theta^{(1)}(t)}$. Thus in contrast to [29] we propose a fit of the full specification of the intensity function including all parameters and the baseline intensity. Furthermore, our aim is to model large networks whereas [29] considered relatively small networks over a long time period $T$. Thus in our asymptotics we let the number of actors converge to infinity instead of $T$. We will argue below (at the end of Section 2.1 and after Assumption (A6) in Section 2.3) that appropriate choices of the censoring factor $C_{n,ij}(t)$ allow for modeling large networks with degrees of the nodes/actors being relatively small compared to the size of the network.

Despite the strong interest in dynamic models for networks, rigorous statistical analysis of corresponding estimators (asymptotic distribution theory) are relatively sparse, in particular in the case of time-varying parameters, as considered here. The temporal models in the literature are usually Markovian in nature. In contrast to that, our continuous-time model based on counting processes allows for non-Markovian structures (i.e. dependence on the infinite past). This increases flexibility in the modeling of the temporal dynamics. Our model also allows for a change of the network size over time without the networks degenerating in the limit. Moreover, we are presenting a rigorous analysis of distributional asymptotic properties of the corresponding maximum likelihood estimators. To the best of our knowledge, no such analysis can be found in the literature, even for the simpler models indicated above.

In Section 2, we discuss our model (Section 2.1), define our likelihood-based estimators, and present our main result on the point-wise asymptotic normality of our estimators in Section 2.2. In Section 3, we demonstrate the finite sample behavior and the flexibility of our approach by presenting an analysis of the Capital Bike-share Data. The proof of our main result is deferred to Section 4. The appendix contains additional simulations where we compare network characteristics as degree distributions, cluster coefficients and diameters of the observed network with networks distributed according to the fitted model. Moreover, we discuss data adaptive bandwidth choices in the appendix.

2 Link-based dynamic models

2.1 Link-based dynamic models with constant parameters

We will first discuss the model described in Section 1.2 with general link function $\Phi$ and *constant *parameter function $\theta\equiv\theta_{0}$. The following form of the log-likelihood for the parameter $\theta$ is shown in [2]:

[TABLE]

where $\Delta N_{n,ij}(t):=N_{n,ij}(t)-N_{n,ij}(t-)$ (and $N_{n,ij}(t-):=\lim_{\delta\to 0,\delta>0}N_{n,ij}(t-\delta)$) is the jump height (either 0 or 1) of $N_{n,ij}$ at $t$. And hence, we obtain the maximum likelihood estimator as

[TABLE]

where $\Theta$ denotes the range in which the true parameter is located. The choice of $\Phi$ as in (1.1) allows for an easy interpretation of the parameters: The intensity has the form $\prod_{k=1}^{q}e^{\theta_{k}X_{n,ij}^{(k)}(t)},$ where $X_{n,ij}^{(k)}(t)$ denotes the $k$-th component of $X_{n,ij}(t)$. Hence, $\theta_{k}$ quantifies the impact of $X_{n,ij}^{(k)}(t)$ on the intensity, given that the remaining covariate vector stays the same.

The presence of the function $C_{n,ij}(t)$ enhances the modeling flexibility significantly. By choosing $C_{n,ij}$ the researcher who applies the model is able to fit the model only to a sub-network. This becomes necessary when it is natural to assume that certain pairs of actors behave fundamentally different from others. For instance, consider a social media network, and contrast pairs impacting each others activities in the network by exchanging messages regularly, with pairs consisting of actors from different social communities hardly interacting with each other. It is intuitive that these two pairs cannot be modeled accurately by the same model. In this situation, it would instead be advantageous to restrict to those pairs who have recently interacted, say, and fit the model only to interactions among them. On the other hand, the interaction intensity of course is dynamic, and thus different pairs might be included over time. This is achieved by the presence of the selector variables $C_{n,ij}(t).$ Also note that $C_{n,ij}(t)=0$ for $t\in[a,b]$ does not necessarily mean that there will be no interactions between $i$ and $j$ in $[a,b]$, it rather means that interactions which happen between $i$ and $j$ in $[a,b]$ are not fitted by our model.

Two things should be noted about these selectors: Firstly, choosing $(C_{n,ij})_{i,j}$ too conservatively is not a problem in the sense that we still estimate the ‘correct’ parameters. For instance, suppose that $C_{n,ij}(t)$ and $\theta$ are the correct quantities to be used in the model (1.1). Assume now that $C^{*}_{n,ij}$ is a predictable selector that is more conservative than $C_{n,ij}$, i.e., $C_{n,ij}(t)=0$ implies $C_{n,ij}^{*}(t)=0$. Then, the observations are given by $t\mapsto N^{*}_{n,ij}(t):=\int_{0}^{t}C^{*}_{n,ij}(s)dN_{n,ij}(s)$ for $t\in[0,T]$. Clearly, $N^{*}_{n,ij}$ is a counting process comprising those jumps of $N_{n,ij}$ at which $C_{n,ij}^{*}$ equals 1. By assumption, $C_{nij}^{*}(t)C_{n,ij}(t)=C_{n,ij}^{*}(t)$ and hence the compensator of $N_{n,ij}^{*}$ is given by $C_{n,ij}^{*}(t)\exp(\theta^{T}(t)X_{n,ij}(t))$. Thus, the processes $N_{n,ij}^{*}$ can also be used to estimate $\theta$. On the other hand, using fewer data of course leads to a loss of information, and this might effect the efficiency of the parameter estimator (cf. Theorem 2.1). We do not attempt here to determine the best $C_{n,ij}$ in a data driven way. Instead, in real data applications, and motivated by this discussion, we attempt to choose $C_{n,ij}$ in a way that is not too liberal. This is illustrated in Section 3, where we set $C_{n,ij}(t)$ equal to zero, if there was no event between $i$ and $j$ for a certain period $\Delta t=(t-\delta,t),$ for some $\delta>0$, so that our model is only fitted to “active” pairs. For pairs with low activity one may look for a different model. Thus a proper choice of $C_{n,ij}(t)$ allows to split up the analysis into different regimes.

Secondly, consider the social media example again. Suppose that the probability that $C_{n,ij}(t)=1$ is the same for all pairs $(i,j)$. We will assume that, with positive probability, links can form at any time $0<t<T$, i.e., $\mathbb{P}(C_{n,ij}(t)=1)>0$. However, it is intuitive that even when more and more people connect to the platform, one particular actor will not acquire an unbounded number of friends. Instead it seems reasonable to assume that its number of friends is bounded. In such a case, the fraction of pairs $(i,j)$ for which $C_{n,ij}(t)=1$ should be of order $\frac{1}{n}$. Hence, it is natural to assume that $\mathbb{P}(C_{n,ij}(t)=1)\to 0$ as $n\to\infty$. This way sparsity in the observations can be captured in the model. For the asymptotic result (Theorem 2.1) to hold, the network cannot be too sparse (essentially an increase in the number of actors must lead to an increase in the number of active pairs). This will be made precise in the assumptions given in Section 2.3.

2.2 Estimation in time-varying coefficient models

In time series applications, it turns out that powerful fits can be achieved by letting the time series parameters depend on time, and this is what we consider here as well. We will use the above model with $\theta$ in (1.1) depending on $t$, or in other words, $\theta=\theta(t)$ is now a parameter function.

An estimator of this parameter function at a given point $t_{0}$ can be obtained by maximizing the following local likelihood function in $\mu$ which is obtained by localizing the likelihood (2.1) for a constant parameter at time $t_{0}$ by means of a kernel $K$

[TABLE]

where $K$ is a kernel function (positive and integrating to one), and $h=h_{n}$ is the bandwidth. The corresponding local MLE is defined as

[TABLE]

with $\Theta$ being the allowed range of the parameter function $\theta$. Recall that we use the following Cox-type form of the intensity:

[TABLE]

With this choice, the local log-likelihood can be written as (up to a term not depending on $\theta$):

[TABLE]

The maximum likelihood estimator $\hat{\theta}_{n}(t_{0})$ studied in this paper is defined as the maximizer of (2.2) over $\theta\in\Theta$, where $\Theta\subseteq\mathbb{R}^{q}$ is an appropriate parameter space. Denote by $L_{n}(t_{0})$ the set of active edges, i.e., the set of all pairs $(i,j)$, such that, $C_{n,ij}(t_{0})=1$. We denote by $|L_{n}(t_{0})|$ the size of the set $L_{n}(t_{0})$. Our main theoretical result, given below, says that for a given $t_{0}\in(0,T)$, the maximum likelihood estimator $\hat{\theta}_{n}(t_{0})$ exists, is asymptotically consistent, and is asymptotically normal.

To formulate our main result, the following technical assumptions are needed.

2.3 Assumptions

Our assumptions do not specify the dynamics of the covariates $X_{n,ij}(t)$ and of the censoring variable $C_{n,ij}(t)$. Instead, we assume that the stochastic behavior of these variables stabilizes for $n\to\infty$. Assumption (A1) is specific to our setting and it states our general understanding of the dynamics, while assumptions (A2), (A3) and (A5) are standard. Assumption (A4) can be found similarly in [29]. It guarantees that the covariates are well behaved. Finally, (A6) and (A7) specifically describe the dependence relations in our context. They quantify the idea that while the network grows the actors get further and further apart and hence influence each other less and less. In what follows, we state our assumptions and briefly discuss their meaning and the intuition behind them.

We denote derivatives by $\partial$. In particular, $\partial_{t}$ and $\partial_{t^{2}}$ refer to the first and second derivative with respect to $t$ respectively. If we derive with respect to a vector $\theta$, $\partial_{\theta}$ refers to the gradient and $\partial_{\theta^{2}}$ refers to the Hessian matrix.

**(A1) Exchangeability **

Assume that for every $n$ and any $s,t\in[t_{0}-h,t_{0}+h],$

the joint distribution of $(C_{n,ij}(t),X_{n,ij}(t))$ is identical for all pairs $(i,j)$, 2. 2.

the conditional distribution of the $q$-dimensional covariate $X_{n,ij}(t)$ given that $C_{n,ij}(s)=1$, has a density $f_{s,t}(y)$ with respect to a measure $\mu$ on $\mathbb{R}^{q}$. This conditional distribution does not depend on $(i,j)$ and $n$. We use the shorthand notation $f_{s}$ for $f_{s,s}$.

The most restrictive part of (A1) is that the conditional distribution of $X_{n,ij}(t),$ given $C_{n,ij}(s)=1,$ does not depend on $i$, $j$. Observe that this holds if the array of $(C_{n,ij},X_{n,ij})_{i,j}$ is jointly exchangeable in $(i,j)$ for any fixed $n$. The additional assumption that the conditional distribution of $X_{n,ij}(t),$ given $C_{n,ij}(s)=1,$ does not change with $n$ is not very restrictive, because it is natural to assume that the distribution depends only on the local structure of the network (Recall the discussion in Section 2.1 in which we assumed that a fixed vertex $i$ has only a bounded number of close interaction partners $j$ while the network grows). We make this additional assumption mainly to avoid stating lengthy technical assumptions allowing to interchange the order of differentiation and integration at several places in the proof.

We add some standard assumptions on the kernel.

**(A2) Kernel and Bandwidth

**Suppose that the Kernel $K$ and the bandwidth $h$ fulfil the following conditions.

$K$* is positive and supported on $[-1,1]$.* 2. 2.

$\int_{-1}^{1}K(u)\mathrm{d}u=1$, $\int_{-1}^{1}K(u)u\mathrm{d}u=0$ and $\max_{-1\leq u\leq 1}K(u)<\infty$. 3. 3.

As $n\to\infty$, $h=o(1)$, $l_{n}:=\frac{n(n-1)}{2}\mathbb{P}(C_{n,12}(t_{0})=1)\to\infty$ with $l_{n}h\to\infty$, and $l_{n}h^{5}=O(1)$.

Note that, $l_{n}$ is the effective sample size at time $t_{0}$, because $\frac{n(n-1)}{2}$ is the number of possible links between vertices, of which, in the average, we observe the fraction $\mathbb{P}(C_{n,12}(t_{0})=1)$. (For directed networks, one simply has to replace $\frac{n(n-1)}{2}$ by $n(n-1)$ in the definition of $l_{n}$.) With this in mind, the assumptions on the bandwidth are standard.

The next assumption states smoothness conditions on the parameter curve $\theta_{0}$.

**(A3) Smoothness of Parameter

**Let $\Theta$ be the convex parameter space and $\theta_{0}:[0,T]\to\Theta$ the parameter function.

$\theta_{0}$* is twice continuously differentiable in a neighborhood of $t_{0}$.* 2. 2.

The value $\theta_{0}(t_{0})$ lies in the interior of $\Theta$.

We continue with some tail conditions on $f_{s}(y)$ and its derivatives. They are fulfilled if, e.g., the covariates are bounded.

(A4) Moment Conditions

For $\mu$-almost all $y$ (where $\mu$ is as in (A1)), $s\mapsto f_{s}(y)$ is twice continuously differentiable. Let $U_{h}:=[t_{0}-h,t_{0}+h]$. There are bounded, open and convex neighborhoods $U$ of $t_{0}$ and $V\subseteq\Theta$ of $\theta_{0}(t_{0})$ such that for all pairs $(i,j)$ and $(k,l)$ and $\tau:=\sup_{\theta\in V}\|\theta\|$,

[TABLE]

For $k\in\{2,3\}$:

[TABLE]

**(A5) Identifiability

$\theta^{T}X_{n,12}(t_{0})=0$*** a.s. (w.r.t. $f_{t_{0}}$) implies that $\theta=0$.

The following assumption addresses the asymptotic behavior of the distributions of the processes $C_{n,ij}(t)$. In particular, for $t$ in a neighborhood of $t_{0}$, we postulate asymptotic stability of the marginal distributions of these processes, and also a certain kind of asymptotic independence of $C_{n,ij}$ and $C_{n,kl}$ for $|\{i,j\}\cap\{k,l\}|=0$.

(A6) Asymptotic Uncorrelatedness I

For $w(u)=K(u)$ and $w(u)=K^{2}(u)/\int K^{2}(v)\mathrm{d}v$ it holds that

[TABLE]

For

[TABLE]

we assume that

[TABLE]

Furthermore, it holds that, as $n\to\infty$

[TABLE]

and, for edges with $|\{i,j\}\cap\{k,l\}|\leq 1$,

[TABLE]

Note firstly that, due to the localization of our likelihood function, time dependence is present only locally around the target time $t_{0}$. Condition (2.10) appears reasonable in our asymptotics where the size of the network increases: Consider, for instance a dynamic social media network, and assume, for example, that we consider data from a certain geographic region. One might assume that at night the number of active pairs, i.e. the pairs with $C_{n,ij}=1,$ is lower than during the day, and that there is a gradual decrease between 8pm and 11pm, say. This time window does not get narrower when $n$ increases and hence a slow change of the distribution over time seems to be a reasonable assumption. Assumption (2.12) holds, for instance, in the following model: Assume that in the previous example communications between pairs end at $\delta_{0}:=8pm$ plus a certain random time $\delta_{n,ij}$, i.e., $C_{n,ij}(t)=\mathbbm{1}(t\leq\delta_{0}+\delta_{n,ij})$. In this case, the ratio of probabilities in (2.12) becomes

[TABLE]

Since we are using a localizing kernel, the length of the interval $[s-\delta_{0},t_{0}-\delta_{0})$ is of the order $h$, and if $\delta_{n,ij}$ has a density, then (2.12) holds.

If we assume that relabeling the vertices does not change the joint distribution of the whole process (i.e. if we assume exchangeability), then, the joint distribution of two pairs $(i,j)$ and $(k,l)$ depends only on $|\{i,j\}\cap\{k,l\}|$. It is thus natural to distinguish the three regimes $|\{i,j\}\cap\{k,l\}|\in\{0,1,2\}$. This pattern will appear again in the next Assumption (A7). Let us for the moment consider $C_{n,ij}$ that are constant over time. Then, in (2.11), the case $|\{i,j\}\cap\{k,l\}|=2$ is true because $\frac{\mathbb{P}(C_{n,ij}=1,C_{n,kl}=1)}{\mathbb{P}(C_{n,ij}=1)^{2}}=\mathbb{P}(C_{n,12}=1)^{-1}=o(n^{2})$ according to Assumption (A2).

We discuss the remaining cases for the uniform configuration model. In this model all vertices have (approximately) the same pre-defined degree $\kappa,$ and we assume that the $C_{n,ij}$ are created as follows: Equip each vertex $i=1,...,n$ with $\kappa$ edge stubs, and create edges by randomly pairing the stubs. After that, discard multiple edges and self-loops. If two vertices $i$ and $j$ are connected after this process, set $C_{n,ij}=1$. We use the same heuristics as e.g. in [28], Chapter 13.2, to compute the probability of edges. Fix $i$ and $j$, then, for any fixed edge stub of $i$, there are $\kappa n-1$ stubs left to pair with, $\kappa$ of which belonging to vertex $j.$ Hence, the probability of connecting to $j$ is given by $\frac{\kappa^{2}}{\kappa n-1}$ as there are $\kappa$ edge stubs from $i$ as well. Thus, for large $n$, we obtain the following probabilities:

[TABLE]

We see now that also for $|\{i,j\}\cap\{k,l\}|\leq 1$, the assumptions (2.11) and (2.13) hold.

The next assumption involves $\theta_{0,n},$ defined as the maximizer of

[TABLE]

where $g$ is defined in (A7). We show later that $\theta_{0,n}$ is uniquely defined, and that $\theta_{0,n}$ is close to $\theta_{0}(t_{0})$ (see Lemma  4.2 and Proposition 4.2, respectively). Define furthermore

[TABLE]

where $X_{n,ij}^{(a)}$ is the $a$-th entry of the vector $X_{n,ij}(s)\in\mathbb{R}^{q}$. Note that, by Assumption (A1), $f_{2}$ and $g$ do not depend on $(i,j)$ and $n$. Keep also in mind that the fact that the covariates $X_{n,ij}(s)$ are vectors implies that $\tau_{n,ij}$, $f_{n,1}$ and $f_{2}$ are matrices and the expectations are to be understood element-wise.

(A7) Asymptotic Uncorrelatedness II

We assume that $f_{n,1}$ depends on $(i,j)$ and $(k,l)$ only through $|\{i,j\}\cap\{k,l\}|$. Moreover, we assume that, for all sequences $(\theta_{n})_{n\in\mathbb{N}}$ with $\theta_{n}\to\theta_{0}(t_{0})$ as $n\to\infty$ and all $u,v\in[-1,1],$ it holds that $f_{n,1}(\theta_{n},t_{0}+uh,t_{0}+vh,(i,j),(k,l))$ converges to a matrix that depends only on $|\{i,j\}\cap\{k,l\}|$. We denote this limit by $f_{1}(\theta_{0}(t_{0}),|\{i,j\}\cap\{k,l\}|)$, and assume that

[TABLE]

For $r_{n,ij}^{(a)}(s),$ we assume that, with $\rho_{n,ijkl}^{(a)}(u,v):=r_{n,ij}^{(a)}(t_{0}+uh)r_{n,kl}^{(a)}(t_{0}+vh)$ and for $|\{i,j\}\cap\{k,l\}|=0$,

[TABLE]

Assumption (A7) specifies in which sense the covariates are asymptotically uncorrelated. For motivating these assumptions build a graph $\mathcal{G}$ with vertices $1,...,n$ and $(i,j)$ being an edge if $C_{n,ij}(t_{0})=1$. Denote by $d_{\mathcal{G}}$ the distance function between edges on $\mathcal{G}$ (that is, the number of edges on a shortest path, i.e., adjacent edges have distance 0). In the same heuristic as given after Assumption (A6), this graph becomes very large (asymptotics over the number of vertices) and sparse ($n$ vertices and of order $n$ edges), because every vertex is incident to at most $\kappa$ edges. In this scenario, the number of pairs of edges $e_{1}$ and $e_{2}$ for which $d_{\mathcal{G}}(e_{1},e_{2})=d$ is of order $(\kappa-1)^{d}\cdot n$, and there are of order $n^{2}$ many pairs of edges in total. Let now $A_{i,j}$ be arbitrary, centered random variables indexed by the edges of $\mathcal{G}$. We make the assumption that $A_{i,j}$ is influenced equally by all $A_{k,l}$ with $(k,l)$ being adjacent to $(i,j)$. In mathematical terms, we formulate this assumption as $\mathbb{E}(A_{i,j}A_{k,l}|d_{\mathcal{G}}((i,j),(k,l))=d)\approx C\cdot\kappa^{-d}$. Then, we obtain for non-adjacent edges $(i,j)$ and $(k,l)$,

[TABLE]

which converges to zero after being multiplied with $l_{n}h\approx nh$ (in this case). Because, in (2.18) and (2.19), we consider only expectations conditional on $C_{n,ij}(t)=1$, we can think of $A_{n,ij}$ being the random variables $\tau_{n,ij}$ (a centered version of it) or $r_{n,ij}^{(a)}$ and the expectations in the above heuristic are conditional expectations conditionally the respective conditions in (2.18) and (2.19). This serves as motivation for these two assumptions. Moreover, unconditionally, $\tau_{n,ij}$ and $\tau_{n,kl}$ (and $r_{n,ij}$ and $r_{n,kl}$) do not need to be uncorrelated.

2.4 The main asymptotic result

Theorem 2.1**.**

Suppose that Assumptions (A1)–(A7) hold for a point $t_{0}\in(0,T)$. Then, with probability tending to one, the derivative of the local log-likelihood function $\ell_{T}(\theta,t_{0})$ has a root $\hat{\theta}_{n}(t_{0})$, satisfying, as $n\to\infty$

[TABLE]

with (denote by $\partial_{t}$ and $\partial_{t^{2}}$ the first and second derivative with respect to $t$ respectively, $\partial_{\theta}$ and $\partial_{\theta^{2}}$ are defined analogously, mind that $\theta$ is a vector)

[TABLE]

and $\tau_{n,ij}(\theta,s)=X_{n,ij}(s)X_{n,ij}(s)^{T}\exp(\theta^{T}X_{n,ij}(s))$ was defined in (2.15). If, in addition, $\frac{|L_{n}(t_{0})|}{l_{n}}\overset{\mathbb{P}}{\to}1$, then $l_{n}$ can be replaced by $|L_{n}(t_{0})|$.

Recalling that (in the case of undirected networks) $l_{n}=\frac{n(n-1)}{2}\mathbb{P}(C_{n,ij}(t_{0})=1)$ is the effective sample size, i.e., the expected number of pairs relevant for estimation of $\theta_{0}(t_{0})$, we see that Theorem 2.1 is a classical asymptotic normality result up to the additional bias term $B_{n},$ which we will discuss next. It holds that

[TABLE]

This is of order O(1) by (2.12) and (2.8). Hence, we get that $B_{n}=O_{P}(1)$. In general, the expectation does not converge to 0. Thus, in general we will have an additional bias term of order $h^{2}$. Let us suppose that one can show $B_{n}-\mathbb{E}(B_{n})=o(1)$ by using some additional assumptions that bound the second moment of this term. We have that $\mathbb{E}(B_{n})=o(1)$ if $\frac{\mathbb{P}(C_{n,12}(t_{0})=0,C_{n,12}(s)=1)}{h\mathbb{P}(C_{n,12}(t_{0})=1)}=o(h)$. This assumption can only hold if only for a negligible minority of edges the membership to the active set changes. In particular, for the extreme case of $C_{n,ij}$ being constant, we have $\gamma_{n,ij}\equiv 0$ and $B_{n}=0$. Hence, the bias term $B_{n}$ is induced by a change in the sparsity of the active set.

Remark 2.2**.**

If one is just interested in consistency, the assumptions can be weakened. In the proof of Theorem 2.1 we need to prove the convergence of a certain quantity to a normal distribution. In order to establish consistency it is sufficient that this quantity converges to zero when being multiplied with a certain null-sequence. In order to show this weaker requirement we do not need the assumptions which impose rates on certain quantities. More precisely we do not need Assumptions (A6), (2.12) and (2.13) and (A7), (2.19). Moreover, the Assumptions (A4), (2.7) and (2.8) may be dropped.

2.5 Direct network modeling

We consider the following general model for the link-based dynamics of a random network, using a multivariate continuous-time counting process approach allowing for arbitrary dependence structure between the links by applying the model for dynamic interactions twice: Once for the formation of new links and once for the deletion of existing links (this separation can also be found in [23]). As before, let $V_{n}=\{1,...,n\}$, be the set of vertices and $L_{n}$ be the set of edges. Note that here we are considering undirected networks. But directed networks can be handled similarly. For a given link $(i,j)$, we let

[TABLE]

Then

[TABLE]

describes the random network, or, equivalently, the (upper half of the) adjacency matrix at time $t$. To describe the dynamics of the links over time we introduce two processes, $N^{+}_{n,ij}(t)$ and $N^{-}_{n,ij}(t),$ counting how often a link $(i,j)$ was added or deleted, respectively, until time $t$. Formally,

[TABLE]

With these definitions, we can write, for $(i,j)\in L_{n}$,

[TABLE]

For $v\in\{+,-\}$, the intensities of the counting processes $N^{v}_{n,ij}(t)$ are here defined as

[TABLE]

with

[TABLE]

for some functions $\gamma^{+}$ and $\gamma^{-}$ respectively, where $\theta^{+}$ and $\theta^{-}$ are two different parameters, determining the addition and the deletion processes, respectively. The vectors $X^{v}_{n,ij}(t)$ for $v\in\{+,-\}$ denote covariates that are assumed to be ${\cal F}_{t}$-predictable. Note that this definition of the intensities makes sure that, as it should be, a link can only be added if it was not present immediately before, and similarly for the removal for a link.

These definitions of the intensities fit into the framework of Section 2.2 with intensity function (2.4), when choosing $\lambda_{n,ij}^{v}(\theta^{v},t)=C^{v}_{n,ij}(t)\cdot\exp(\theta^{v}(t)^{T}X^{v}_{n,ij}(t))$ with $C_{n,ij}^{v}(t)$ being predictable $\{0,1\}$-valued processes that fulfill $C_{n,ij}^{+}(t)=0$ if $Z_{n,ij}(t-)=1$, and $C_{n,ij}^{-}(t)=0$ if $Z_{n,ij}(t-)=0$. Again, as in Section 2.2, we allow that the parameter is a function of time. To sum it up: The processes $N_{n,ij}^{+}$ are modeled with intensity $\lambda_{n,ij}^{+}(\theta_{0}^{+},t)$ and the processes $N_{n,ij}^{-}$ are modeled with intensity function $\lambda_{n,ij}^{-}(\theta_{0}^{-},t)$. Our model allows the covariates $X_{n,ij}^{v}$ and the true parameter functions $\theta_{0}^{v}$ to be different for $v=^{\prime}+^{\prime}$ and $v=^{\prime}-^{\prime}$. For estimating the parameters, we consider observations of the same type only, i.e., we will compute two maximum likelihood estimators: the estimator of $\theta_{0}^{+}(t)$ based on the processes $N_{n,ij}^{+}$, and the estimator for $\theta_{0}^{-}$ based on the processes $N_{n,ij}^{-}$. Both estimators can be treated as coming from an interaction based model and hence the theory from Section 2.2 can be applied.

3 Application to Bike Data

Here we illustrate the finite sample performance of our estimation procedure described above, by considering the Capital Bikeshare (CB) Performance Data, publicly available at http://www.capitalbikeshare.com/system-data. The available data describes the usage of the CB-system at Washington D.C. from 2010 to 2018. However, for computational reasons, and in order to keep the presentation concise, we will present two analyses of sub data sets. In the first analysis, we study the bike data from Jan 2012 to March 2016. In order to reduce computational complexity we aggregate the data over days. In this first analysis it is our main interest to predict the activity of an edge based on the past. In the second analysis we focus on a short period, April and May 2018, and we keep the time-continuous scale of bike events.

The available data set does not contain bike rentals over several days or below 60 seconds, and service rides are excluded. While the last two seem not restrictive, for bike rentals for more than a day, we note the pricing structure of CB: No matter which pass you buy (single ride, day pass, 3-day pass, 30-day pass or annual membership) the basic fee always includes only bike rides for less than 30min. If you keep a bike for longer than 30min extra fees apply (if you buy e.g. a day pass for $8 you have an unlimited number of bike rides up to 30min, but if you keep the bike for 10h, you will have to pay$142). We assume therefore that such long rentals are not the companies primary business and that they occur only very rarely. It is however possible to return a bike to a bike station and immediately re-rent it. Such that, in practice, if you want to make a way on the bike which takes you more than 30min, you can make an intermediate stop to avoid cost. Lastly, we have no information about the status of the stations themselves. In particular, we do not know if a station is empty or full.

It should be noted that, while we believe that this example serves as a serious and interesting illustration of our proposed method, it is not meant to be a full-fledged analysis of bike sharing performance. We would rather like to make the case for the potential of the model along with the estimation strategy presented in this paper by arguing that intuitively convincing results for the bike sharing data indicate that the model might also be beneficially used in more complex situations (i.e., without a strong a priori intuition).

Generally, in both analyses we consider the bike stations to be vertices in the network. Whenever somebody rents a bike at station $i$ and returns it at station $j$, we consider this an event from $i$ to $j$ and in this case we say that $(i,j)$ has been used. In May 2018 the CB network comprises 527 bike stations which were used at least once in April or May 2018. This results in a total of 277,202 possible directed combinations. Of these 277,202 directed combinations, only 39,722 connections have been used at least once in April 2018, and only 9,131 combinations have been used ten times or more. We conclude that the network is very sparse and that it is very challenging for a model to capture the entire biking behavior among all 527 bike stations. Thus, we restrict our analysis to some subset of pairs of bike station which (we assume) can be reasonably modeled by the same model. Consider an example: Bike stations in Alexandria and Derwood are 50km apart, and, on the other hand, some bike stations in downtown Washington are just separated by one block. Certainly, every now and then, somebody might take such a bike ride, but we cannot expect that our model will capture all these special cases. This restriction is realized by appropriate choice of the indicator functions $C_{n,ij}$ (see below for more details). Note here that, on a general level, $C_{n,ij}(t)=1$ means that the pair $(i,j)$ is, at time $t$, regarded as being part of the model in the sense that events from $i$ to $j$ can be captured by our model. On the other hand $C_{n,ij}(t)=0$ does simply mean that the pair $(i,j)$ does not belong to those edges of interest to us. While there might still be bike rides from $i$ to $j$, we do not attempt to model them by using our model.

3.1 Analysis 1: January 2012 till March 2016

In this part of the analysis we intend to model the biking activity on one weekday (Friday) based on the past. By biking activity on an edge $(i,j)$ we mean the number of bike rides between bike stations $i$ and $j$. Direction does not matter to us in this part of the study. The decision to only model one weekday was made mainly to reduce the computational burden. Comparing the results for different weekdays might be instructive, in particular comparing a regular working day and a day on the weekend. Note also that in the period of four and a quarter years, which we consider, eight Fridays were actually public holidays (thus being possibly more like a weekend than a weekday). They were Independence Day (2013, 2014), the Friday after Thanksgiving (2012, 2013) and Fridays during Christmas and New Year’s Holidays. Notice, however, that the parameter function is allowed to change over time. Thus, we assume that the influence of public holidays is not causing any problems.

Figure 1 shows some summary statistics of the data. In Figure 1(a), we see the number of available bike stations, which is strongly increasing. Figure 1(b) shows the number of bike tours on Fridays. An obvious seasonal periodicity is visible with low activity in winter. In order to reduce computational complexity to a minimum (fitting the model takes several minutes on standard laptop), we assume that the covariates change only at midnight and stay constant over the day. Furthermore, we estimate the time-varying parameter function $\theta$ only for one time point per day, namely 12pm noon. The next paragraph contains more details.

Since we do not consider any asymptotics here, we omit the index $n$. Time $t$ is measured in hours of consecutive Fridays. So, if $k$ is the current week, and $r$ is the time on Friday (in 24h), then $t:=(k-1)\cdot 24+r$. Thus, with $r_{t}:=(t\mod 24),$ the quantity $k_{t}:=\frac{t-r_{t}}{24}+1$ gives the week the time point $t$ falls into. The processes $N_{i,j}(t)$, counting the number of tours between $i$ and $j$ on Fridays, are modeled as counting processes with intensities $\lambda_{i,j}(\theta(t),t):=\alpha(t)\exp(\theta^{T}X_{i,j}(k_{t}))\cdot C_{i,j}(k_{t})$. The covariate vector $X_{i,j}(k_{t})$ and the censoring indicator $C_{i,j}(k_{t})$ will be defined later. Note that they both only depend on $k_{t}$, i.e. on the current week, and not on the actual time on the Friday under consideration. The function $\alpha$ is 24 periodic and integrates to one over a period, i.e., $\alpha(t)=\alpha(t+24)$ and $\int_{t}^{t+24}\alpha(s)ds=1$. The role of the (unobservable) function $\alpha$ is to argue that discretizing the biking activity is not introducing a bias even when the biking activity varies over the day. Suppose now, that our target is the estimation of the parameter vector $\theta(t_{0})$ with $t_{0}=(k_{t_{0}}-1)24+r_{0}$ and $r_{0}=12$, say. We choose a piecewise constant kernel $K$ with $K((24k+x)/h)=K(24k/h),$ for all $k\in\mathbb{N}$ and $0\leq x<24$. Substituting in these choices of the intensity and the kernel to the log-likelihood (2.2), we see that our maximum likelihood estimator maximizes the function

[TABLE]

where $\int_{k\cdot 24}^{(k+1)\cdot 24}dN_{i,j}(t)$ gives the number of tours between $i$ and $j$ on the Friday in week $k,$ and where $K_{\kappa}(k)=K(k/\kappa)$ with $\kappa=h/24$. In our empirical analysis, we chose $K_{\kappa}(k)$ as triangle weights with support $\{-\kappa,...,\kappa\}$ and considered only integer choices of the bandwidth $\kappa$. The bandwidth choice is discussed at the end of this section. Note that due to this discretization we essentially obtain a sequence of generalized linear Poisson models with time varying parameters. In the second analysis, in Section 3.2, we use the full time-continuous potential of the model for dynamic interaction networks.

We explain now the choice of our covariate vector $X_{i,j}$. Denote by $\Delta_{i,j}(k,d)$ the number of tours between $i$ and $j$ on day $d$ in week $k$, where $d=4$ means Monday and $d=7$ refers to Thursday (for us the week starts on Fridays, i.e. Friday is $d=1$). For $r\in(0,1)$, we encode the activity between $i$ and $j$ in week $k$ as $A_{i,j,k}=(1-r)\sum_{d=4}^{7}r^{7-d}\Delta_{i,j}(k,d)$ (mind the limits of the summation - Fridays are not included). In our simulations, we chose $r=0.8$ (this choice is somewhat arbitrary, and a full study of the data would include investigating the sensitivity of the parameter estimate on the choice of $r$ as well as a data driven choice. We do not attempt to do this here). We construct a network $G(k),$ for every week $k,$ by connecting $i$ and $j$, if and only if, there was at least one tour on the Friday in that week. We denote by $I_{i,j,k}$ the number of common neighbors of $i$ and $j$ in the graph $G(k)$. We let $d_{i,k}$ be the degree of node $i$ in $G(k)$, $T_{i,j,k}$ the number of tours between $i$ and $j$ on the Friday in the $k$-th week, and $T_{i,j,k,k-1}=(T_{i,j,k}+T_{i,j,k-1})/2$ the average number of tours on the two Fridays in weeks $k$ and $k-1$. Finally we collect everything in the covariate vector:

[TABLE]

The censoring indicator function $C_{i,j}$ is defined to be equal to zero, if there was no tour between stations $i$ and $j$ in the last four weeks. In other words, we attempt to model only those connections which are used regularly in the considered time frame. In summary, we estimate a total of six parameter curves, corresponding to the effects of six covariates in our model:

• •

$\theta_{1}(t)\quad\triangleq\quad$ baseline

• •

$\theta_{2}(t)\quad\triangleq\quad$ activity between stations on previous week-days

• •

$\theta_{3}(t)\quad\triangleq\quad$ common neighbors of stations

• •

$\theta_{4}(t)\quad\triangleq\quad$ popularity of station, measured by degrees

• •

$\theta_{5}(t)\quad\triangleq\quad$ activity between stations on two previous Fridays

• •

$\theta_{6}(t)\quad\triangleq\quad$ inactivity between stations on two previous Fridays

Figure 2 gives some impression of the distribution of the covariates over time. We only consider covariates 2-5 (that is the entries 2-5 of the covariate vector $X_{i,j}(k)$). The first covariate is always equal to one, while the last one is an indicator and thus either zero or one, and so we do not give plots for them. Each panel in Figure 2 shows the 50%, 80%, 90% and 99% quantiles of the respective covariate. We see that the quantiles mainly stay on a moderate level with some larger values in between. This effect is more pronounced for the activity based covariates 2 and 5.

The resulting estimated parameter curves are shown in Figures 3 and 4. All calculations have been executed on the BwForCluster (cf. Acknowledgement). Since we expect the parameter function to vary slowly, we used the last estimated value as initial value for the estimation at the next point in time. In all six parameter curves in Figures 3 and 4, the solid curves show the estimated parameter curve and the dotted curves indicate approximative 95% point-wise confidence sets, which we obtained by omitting the bias in Theorem 2.1 and approximating $\Sigma$ at time $t_{0}$ by $\frac{1}{|L_{n}(t_{0})|}\partial_{\theta^{2}}\ell_{T}(\hat{\theta}_{n}(t_{0}),t_{0})$, where $|L_{n}(t_{0})|$ is the number of active edges at time $t_{0}$. In all plots we observe a clearly visible seasonality. Looking at Figure 3(b), we see that activity during the week (Monday to Thursday) is more important during the winter months than in the summer. A plausible interpretation for this might be that the opportunist cyclists might be less active in winter because of the colder weather. So only those keep using a bike, who ride the same tour every day regardless of the weather. This makes the activity during the week a better predictor.

Figure 3(c) shows that the number of common neighbors always has a significant positive effect on the hazard. This reflects the empirical finding that observed networks cluster more than totally random networks (e.g. see [15]).

The influence of the popularity of the involved bike stations is investigated in Figure 4(a) (measured by the degree of the bike station). Interestingly, it always has a significant negative impact. The size of the impact is higher in the summer months, which again supports the hypothesis that in summer the behavior of the network as a whole appears more random than in winter. But still, the negative impact is a bit unforeseen. This finding can be interpreted as the observed network having no hubs. Another reason for this effect might be, that stations can only host a fixed number of bikes: If a station $i$ is empty, no new neighbors can be formed. A similar saturation effect happens if a lot of bikes arrive at station $i$. Moreover, it is plausible that effects caused by the degrees are already included in 3(b), as well as in Figure 4(b). They show the effect of the bike rides on the days immediately preceding the current Friday, and the effect of the average number of bike tours on the last two Fridays, respectively. In Figure 4(b), we observe a similar behavior as in Figure 3(b) (even more pronounced): In summer the predictive power of the tours on the last two Fridays is significantly lower than in winter, underpinning the theory that the destinations in summer tend to be based on more spontaneous decisions. Finally, in Figure 4(c), we observe that no bike tours on the last two Fridays between a given pair of stations always has a significant negative impact on the hazard. Again a very plausible finding.

We are currently working on testing whether the parameter functions depend on time, i.e., on testing for constancy of the parameter functions. For a complete data analysis it would then be interesting to add time as a covariate (or time dependent covariates), and to see if the parameter functions show always a significant time-dependency.

Modeling other network characteristics. In stochastic network analysis, a central strand of research is concerned with the question of whether characteristics observed in real networks can be adequately mimicked by stochastic network models. Important characteristics are degree distribution, clustering coefficient and diameter (these and other characteristics can be found in [15] Chapter 2.2, we define them also in the appendix). As in [41], Chapter 4, we compare these three characteristics with a typical network produced by our model. In order to see how much our fitted model is able to capture these characteristics, we have simulated 3840111We chose to simulate 3840 networks, because we had $32$ cores available, and on each of the cores we ran 120 predictions, which could be done in reasonable time. networks corresponding to three randomly chosen days, by using the network model with the fitted parameters of the corresponding day. We then compared the simulated three characteristics on these three days to the ones observed in the networks (this way of assessing the goodness of fit is also used in [13]). Here, we present the results for the degree distribution on 7th December 2012. The other results are reported in the appendix.

In our analysis, we consider fitting sub-networks defined by the popularity of their edges: For given values $0\leq l_{1}<l_{2}\leq\infty$, the network is constructed by placing an edge between a pair of nodes $(i,j),$ if the number of tours between $i$ and $j$ falls between $l_{1}$ and $l_{2}$. Different ranges of $l_{1}$ and $l_{2}$ are considered. The idea is to consider the network of low frequented tours (for $l_{1}=1$ and $l_{2}=3$) up to the network of highly frequented tours (for $l_{1}=10$ and $l_{2}=\infty$).

Figure 5 shows the simulated degree distributions for six different choices of $l_{1}$ and $l_{2}$. The dotted lines indicate 10% and 90% quantiles of the simulated graphs, and the solid line shows the true degree distribution. We see that, in all six cases, the approximation is reasonable accurate, in particular if one takes into account that we did not specifically aim at reproducing the degree distributions. The plots show that the largest degree of the simulated networks and the observed network lie not too far from each other, and the overall shape of the degree distribution is captured well. It should also be noted that we used only six covariates, whereas in other related empirical work much higher dimensional models have been used, see e.g. the discussions in [29].

Brief remark on choice of bandwidth via one-side cross validation. To choose the bandwidth, we calculate a local linear estimate with a one-sided kernel $K_{+,\kappa}(k)=K_{\kappa}(k)\mathbbm{1}(k<0)$. For all values of $\kappa$, the fitted value of the conditional expectation of $X_{i,j}(k_{t_{0}}),$ given the past, is compared with the outcome of $X_{i,j}(k_{t_{0}})$. This is done for all non-censored edges. The results for different bandwidths are shown in Figure 6. We see that the prediction error of the model decreases, until we reach the bandwidth $\kappa=23$. In one-sided cross-validation, one now makes use of the fact that the ratio of asymptotically optimal bandwidths of two kernel estimators with different kernels, $K$ and $L$ is equal to $\rho=[\int K^{2}(u)\mathrm{d}u(\int u^{2}L(u)\mathrm{d}u)^{2}(\int L^{2}(u)\mathrm{d}u)^{-1}(\int u^{2}K(u)\mathrm{d}u)^{-2}]^{1/5}$. For a triangular kernel, and its one-sided version, we get $\rho\approx 1.82$. The one-sided CV bandwidths is given by dividing 23 by $\rho$ which yields bandwidth roughly twelve (here we also only consider integer bandwidths). More details on the one-sided cross-validation approach are presented in the appendix.

3.2 Analysis 2: May 2018

In this analysis, we study the biking behavior in May 2018 in more detail. In particular, we can (in contrast to before) assume that the number of bike stations remains constant over the observation period. Our main interest in this part lies in illustrating how the model can be used to understand how the system would change if another bike station were built. Let us firstly look at the distribution of bike rides over four weeks in April: this is shown in Figure 7. We see a clear daily pattern: During weekdays the number of bike rides spikes in the morning and in the afternoon while it shows deep valleys (going almost down to zero) at night and not so deep valleys around midday. The weekends show a clearly different pattern by not exhibiting the morning/afternoon spikes so visibly. The only weekdays which depart from these pattern are April 24 and 25. These were both rainy days (we use weather data from the weather station at Washington D.C. Dulles Airport, as reported on Weather Underground). However, we should say that there are other rainy days which do not show such a visible effect. Interestingly, April 16th (Emancipation day that year, a public holiday) is showing similar behavior as the other weekdays but with a smaller number of bike rides (maybe one would have rather expected that public holidays behave like weekends).

In this analysis we choose to restrict to those 9,131 connections which have been used ten times or more in April. However, for this part, we keep the directions of the bike rides. Thus, for any two bike stations, in our model, we let $C_{n,ij}$ be the indicator that the directed connection $(i,j)$ has had ten or more rides in April. The covariates which we use here are based on the distances between two bike stations and their densities. The distance from station $i$ to $j$ is given by the time it takes to go from station $i$ to station $j$ on a bike. These times were computed by using Google maps on a weekday afternoon. Note that the travel time from $i$ to $j$ can be different than that from $j$ to $i$ because of different one-way structures of the streets and possible ascends. Sometimes people on bikes do not quite follow traffic regulations but nevertheless we assume that the travel time is a reasonable measure of distance. We will denote the distance from station $i$ to station $j$ by $d_{i,j}$. The density of a bike station is measured in terms of the number of neighboring stations. We denote by $n(i)$ the average of: 1) the number of bike stations which can be reached from $i$ in less than three minutes and 2) the number of bike stations from which $i$ can be reached in less than three minutes. Our intuition was that if a bike station is full or empty, then people would have to go to another bike station instead. However, we assume that people would not accept an arbitrarily long detour. Therefore, we chose the limit of three minutes bike riding time (not walking time) for neighboring stations. Of course, this threshold is somewhat arbitrary and a full analysis should consider the sensitivity of the results with respect to this threshold. Note lastly, that people can see the availability of bikes and empty docks in advance. Thus, in our model we assume that the bias introduced by people arriving at a full station and being forced to go to another station (no matter how far away) is negligible.

Lastly, we mention that we did not include the precipitation as a covariate for two reasons: Firstly, we wanted to have an hourly analysis but we only had daily data for precipitation. Thus, we could not determine the times of actual rain and, probably, prospective rain in the evening is not going to impact the biking activity in the morning: In a CB member survey from 2016 a bit more than half of the respondents said that one of their main reasons for joining CB is to have access to one-way trips, thus we assume that the possibility of rain in the evening would not stop people from using a bike in the morning. Secondly, and possibly more severely, it is difficult to include covariates which are constant across all connections. If such a covariate were zero it would mean that its corresponding parameter has no influence at all on the intensities. While in theory identification is possibly still valid, the practical computation will break down.

Our aim is to use the model in order to quantify the possible impact of a new bike station on the system. We let $h\approx 1.1$hours (this bandwidth was chosen by the same procedure as outlined at the end of the precious subsection). With such a short bandwidth we will not smooth out differences between morning and afternoon. The covariate vector $X_{n,ij}$ is given by

[TABLE]

Note that, in order to avoid taking the logarithm of zero, all quantities have been bounded from below by 1. In Figure 8 we show the estimated parameter values for the second week of May. The solid lines show the estimates while the dotted lines show approximative 99% point-wise asymptotic confidence regions as provided in the theory above (they were approximated in the same way as in the previous section by the Hessian of the likelihood at the estimate). We assume that the bias is negligible. The results for the other weeks look similarly. Therefore, we consider the results for the entire month only for the intercept and the two covariates indicating the number of neighbors of starting and ending station, cf. Figures 9-11.

We summarize some observations about these estimates:

The estimates of the intercept shows valleys over night and summits during the day. On May 16th-18th and 22nd the intercept is lower than usual. Most weekdays show a valley around midday between two peeks in the morning and in the afternoon (very pronounced on Friday, 18th, and Wednesday 23rd to Friday 25th). Monday, 28th, doesn’t show this behavior. On weekends these peeks are sometimes visible, sometimes not. 2. 2.

On weekdays: Estimates of the parameter corresponding to the log-number of bike stations in the neighborhood around the ending station show negative valleys towards the afternoon/evening and become positive during night time. 3. 3.

On weekdays the log-number of neighboring stations around the start station show the opposite behavior: Estimates show negative valleys towards the morning. Sometimes they become positive during night time.

Before discussing the interpretation of these findings we mention findings about the second and third covariate: The first three covariates model the influence of the log-time ($\log(d_{i,j})$) on the intensity of bike rides as a quadratic function when the densities of start and end station remain fixed. The estimates from Figure 8 suggest that the third covariate, that is the factor in front of $\log(d_{i,j})^{2}$, is often negative. That means that the parabola is open to the bottom and has hence a maximum. As the parameters change over time, the location of this extremum is also moving. In Figure 12 the location of the extremum of this parabola is shown. We observe the following:

The locations of the extrema lie almost always well in the region between 1 and 5 min. 2. 5.

There is no clear pattern visible. 3. 6.

The location of the maximum is not changing much.

A possible interpretation of these findings is as follows: Point one is a very plausible observation as people use bikes less during night. The exceptions mentioned in point one are rainy days (cf. Figure 13). It is interesting to note that other rainy days (like May 6th) are not so much visible in the estimates. We should note that we only have average precipitation information available per day. So it might be possible that the rain came during night when there was not much biking anyway. The peeks in the morning and in the afternoon correspond to the increased activity in the morning and in the afternoon which we observed earlier. Note that May 28th was Memorial Day in 2018. Thus it is not surprising that this day does not show these peeks which might correspond to people commuting to work. The second and third observations mean that in the evening a larger number of bike stations close to the destination yields a lower intensity while in the morning a larger number of bike stations close to the origin yields a lower intensity (however, this effect is not so pronounced). We explain this observation by using commuters: Commuters start their way to work in Washington D.C. possibly from central locations (like bike stations close to Union Station or Metro Stations) and disperse from there through the city to their respective work places. Thus, such stations might empty in the morning and people have to walk to other bike stations. This would yield an effect as observed: Stations with many near by stations in the morning share their traffic with the remaining stations. In the evening the reverse effect is happening: People return from their non-central work places to central bike stations causing these bike stations to fill up and hence forcing people to go to other empty bike stations near by.

During night time we observed a positive effect of the number of neighbors of the destination station. A possible explanation for this might be that the stations were most likely built in a way such that the density of stations is large if the demand for bikes in this neighborhood is large. Thus, the number of neighboring bike stations serves as a proxy for popularity of a bike station. Keeping this in mind, we interpret the positive parameter during night time as indicator for the hypothesis that, at night, people prefer to go to central locations (like metro stations) by bike. The absence of positivity of the number of neighbors of the origin station indicates that the reverse is not always happening: People do not leave from central locations (this is plausible as we would expect more people to go home at night).

The findings four to six are mostly interesting because they indicate that an actual parabola is fit indicating a strict convex behavior of the intensity. It is important to note that by choosing the parameters two and three the model fit could result in a strictly monotone fit of the parabola (when the extremum of the parabola is located outside the interval $[0\textrm{min},30\textrm{min}]$). That this is not the case indicates that there is a non-proportional change in the activity as the distance between bike stations changes. Generally, the location of the maximum in the short distance range indicates that people prefer to take shorter routes rather than longer routes (cf. also the discussion of the CB pricing system).

If CB was to built a new bike station we would assume that they do it because in that region they suspect high biking activity (e.g. because bike stations are constantly empty or full). In that sense we assume that the idea that a high density of biking stations indicates a region interesting for biking is not violated. However, adding a new station changes the covariates in the network and thus we can use the parameter estimates from above in order to predict how the biking activity at other bike stations might change. The introduction of a new station changes for example the number of neighboring bike stations of its neighbors. Thus we could predict how much traffic is diverted from the existing stations to a new station. Moreover, we can also use it to predict the number of bike rides to the station and from the station. This can help to find an accurate size for the new bike station.

In this example we saw that a time varying parameter choice is useful in order to be able to distinguish morning and afternoon as well as rainy and non-rainy periods. Moreover, we illustrated how the covariates could be chosen in order to assess the effect of adding a new bike station.

4 Proof of Theorem 2.1

In the proof, we do not distinguish explicitly directed and undirected networks: in the undirected case, we always assume $i<j$, moreover we will need $l_{n}=O(n^{2}\mathbb{P}(C_{n,12}(t_{0})=1)$, which is true in both cases. The processes $N_{n,ij}$ are counting processes with intensity given by $\lambda_{n,ij}(\theta_{0}(t),t)$. We can decompose these counting processes as (Doob-Meyer Decomposition, see e.g. [2] Chapter II.4)

[TABLE]

where $M_{n,ij}$ is a local, square integrable martingale. We use this decomposition of the counting processes in order to decompose the likelihood and its derivatives. Let $P_{n}(\theta)$ be defined as

[TABLE]

Note that we do not make the dependence of $P_{n}(\theta)$ on $t_{0}$ explicit in the notation. In order to reduce notation, we write for the derivative of a function $\psi(\theta)$ of one variable $\theta$ (which might be a vector) simply $\psi^{\prime}:=\partial_{\theta}\psi$ and $\psi^{\prime\prime}:=\partial_{\theta^{2}}\psi$. Using $P_{n}(\theta),$ we can write

[TABLE]

Recall that $\theta_{0,n}$ is defined as the maximizer of $\theta\mapsto\int_{0}^{T}\frac{1}{h}K\left(\frac{s-t_{0}}{h}\right)g(\theta,s)ds$, where $g$ is defined in (A7). Note that the function $g$ does not depend on $n$, see Assumption (A1). Lemma 4.2 shows that $\theta_{0,n}$ is uniquely defined. The value $\theta_{0,n}$ is the deterministic counterpart of the random quantity $\widetilde{\theta}_{n}(t_{0})$ that is defined as the solution of $P^{\prime}_{n}(\widetilde{\theta}_{n}(t_{0}))=0$. The existence of the latter is considered in Proposition 4.3.

In all lemmas and propositions of this section, we assume that Assumptions (A1), (A2) and (A5) hold as they are permanently used (also implicitly in other assumptions). The other assumptions will be mentioned in those places where they are needed.

Lemma 4.1**.**

We have

[TABLE]

Equality holds, if and only if, $\theta_{0}(s)^{T}y=\theta^{T}y$. In particular, $\theta_{0}(s)$ is the unique maximizer of $\theta\mapsto g(\theta,s)$.

Proof.

Note that, for arbitrary $y\in\mathbb{R},$

[TABLE]

implies that the differentiable function $x\mapsto xe^{y}-e^{x}$ has the unique maximizer $x=y$. This also implies the second statement of the lemma by (A5). ∎

Fact 4.1**.**

Assume (A4), (2.6) and (A3) hold. For $j\in\{0,1,2\},k\in\{0,1,2,3\},$ with $j+k\leq 3,$ the partial derivatives of order $j$ of the function $g(\theta,s)$ with respect to $s,$ and of order $k$ with respect to $\theta,$ exist, for $(t,\theta)\in U\times V$ (cf. Assumption (A4) for a definition of $U$ and $V$). The partial derivatives can be calculated by interchanging the order of integration and differentiation in (2.17). All these partial derivatives of $g(\theta,s)$ are absolutely bounded on $U\times V$. For the calculation of the first two derivatives of $g$ with respect to $\theta,$ differentiation and application of the expectation operator can be interchanged in (2.3). The matrix $\Sigma$ is invertible.

Proof.

The statement of this fact follows immediately from (2.6) of Condition (A4). Note that the functions $\theta_{0}$, $\theta_{0}^{\prime}$ and $\theta_{0}^{\prime\prime}$ are absolutely bounded in a neighborhood of $t_{0}$. This holds because these functions are continuous in a neighborhood of $t_{0}$, see (A3). Invertibility of $\Sigma$ is a consequence of (A5). ∎

Lemma 4.2**.**

Assume Fact 4.1 holds and that $\Theta$ is convex. For $n$ large enough, $\theta_{0,n}$ (the maximizer of (2.14)) is well defined and unique. It holds that $\theta_{0,n}\to\theta_{0}(t_{0})$ as $n\to\infty$. In particular, $\theta_{0,n}\in V$ for $n$ large enough.

Proof of Lemma 4.2.

The function $\theta\mapsto g(\theta,t_{0})$ is strictly concave and hence $\theta_{0}(t_{0})$ is its unique maximizer (cf. Lemma 4.1). Thus

[TABLE]

is strictly concave too. Moreover, we know that $\partial_{t}g(\theta,t)$ is absolutely bounded on $U\times V$. This implies that $\phi_{n}$ converges to $g(\theta,t_{0})$, uniformly on $V$. Hence, $\phi_{n}$ has a local maximizer $\theta_{0,n}$ in the open set $V$. By strict convexity, $\theta_{0,n}$ is the unique global maximum. The convergence $\theta_{0,n}$ to $\theta_{0}(t_{0})$ follows by uniform convergence of $\phi_{n}$ to $g$. ∎

Lemma 4.3**.**

Assume Fact 4.1 holds. With $\Sigma_{n}=-\int_{-1}^{1}K(u)\int_{0}^{1}\partial_{\theta^{2}}g(\theta_{0}(t_{0})+\alpha(\theta_{0,n}-\theta_{0}(t_{0})),t_{0}+uh)\mathrm{d}\alpha\mathrm{d}u$, we have

[TABLE]

Moreover, the sequence

[TABLE]

is bounded, and it holds that $v_{n}\to v$, as $n\to\infty$.

Proof.

Using Lemmas 4.2 and Fact 4.1, we conclude that the integrand

[TABLE]

(note that $u\in[-1,1]$ and $\alpha\in[0,1]$). The first statement of the lemma follows by an application of Lebesgue’s Dominated Convergence Theorem, and the fact that $\partial_{\theta^{2}}g$ is bounded as a continuous function on a compact set. The second statement of the lemma follows similarly. ∎

Proposition 4.2**.**

Assume Fact 4.1 holds. We have, for $t_{0}\in(0,T),$

[TABLE]

Proof of Proposition 4.2.

Since $\theta_{0}(s)$ maximizes $\theta\mapsto g(\theta,s)$ (cf. Lemma 4.1), we have $\partial_{\theta}g(\theta_{0}(s),s)=0$. Furthermore, by definition of $\theta_{0,n},$ we have

[TABLE]

Having observed that, we compute, for $h$ small enough,

[TABLE]

$\Sigma_{n}$ converges to the invertible matrix $\Sigma$ by Lemma 4.3. The first integral is of order $h^{2}$. This follows by a Taylor expansion in the time parameter:

[TABLE]

By Lemma 4.3, $v_{n}$ is bounded. Thus, together with (4.6), we have established

[TABLE]

The statement of the proposition now follows from $v_{n}\to v$. ∎

Lemma 4.4**.**

Assume Fact 4.1, (A4) (2.6), (2.9), (A6) (2.10), (2.11) and (A7), (2.18) hold. We have

[TABLE]

For any $k,l\in\{1,...,q\},$ it holds that

[TABLE]

Moreover,

[TABLE]

where $P_{n}^{{}^{\prime}(r)}$ denotes the $r$-th component of $P^{\prime}_{n}$, the supremum runs over $k,l,r\in\{1,...,q\},$ and $\theta\in V$.

Proof.

We start by showing that $P_{n}^{\prime}(\theta_{0,n})=o_{P}(1)$. This holds, if $\mathbb{E}(\|P_{n}^{\prime}(\theta_{0,n})\|)=o(1)$. Define $\rho_{n,ij}(\theta,s):=\|X_{n,ij}(s)\|\cdot\left|\exp(\theta_{0}(s)^{T}X_{n,ij}(s))-\exp(\theta^{T}X_{n,ij}(s))\right|$. By positivity of $\rho_{n,ij}(\theta,s),$ we may apply Fubini’s Theorem, and thus we compute

[TABLE]

The expectation in the integral expression can be bounded by applying a Taylor expansion:

[TABLE]

where $s_{u}=t_{0}+uh$. Now, by (2.9) in Assumption (A4), the expectation in the last upper bound is bounded by a constant $C,$ uniformly in $u\in[-1,1]$. Using $\sup_{u\in[-1,1]}\|\theta_{0}(t_{0}+uh)-\theta_{n,0}\|=o(1),$ we obtain

[TABLE]

where the last equality is a consequence of (2.10). This shows (4.7).

We now show (4.8). With $\partial_{\theta^{2}}g(\theta_{0,n},s)=-\mathbb{E}(\tau_{n,ij}(\theta_{0,n},s)|C_{n,ij}(s)=1),$ Fact 4.1 gives

[TABLE]

For (4.8), it suffices to show:

[TABLE]

For the proof of (4.11), we note that with $a_{n}(u)=\frac{\mathbb{P}(C_{n,12}(t_{0}+uh)=1)}{\mathbb{P}(C_{n,12}(t_{0})=1)},$

[TABLE]

Here we use (2.10), and $\theta_{0,n}-\theta_{0}(t_{0})=o(1)$ (see Proposition 4.2).

For the proof of (4.10), we write $K_{h,t_{0}}(s):=K\left(\frac{s-t_{0}}{h}\right)$ and

[TABLE]

We will apply Markov’s inequality to show that this term converges to zero. When squaring the above sum, we can split the resulting double sum into three parts, depending on whether $|\{i,j\}\cap\{k,l\}|=0,1$ or $2$. Thus we have to show that the following three sequences converge to zero:

[TABLE]

where $\kappa_{n,ij}(\theta_{0,n},s):=-C_{n,ij}(s)\tau_{n,ij}(\theta_{0,n},s)+\mathbb{P}(C_{n,ij}(s)=1)\partial_{\theta^{2}}g(\theta_{0,n},s),$ and $\bar{\kappa}_{n,ij}(\theta_{0,n}):=\int_{0}^{T}K\left(\frac{s-t_{0}}{h}\right)\kappa_{n,ij}(\theta_{0,n},s)ds$. Now note that

[TABLE]

and that the sum in (4.12) has $O(n^{2})$ terms, (4.13) comprises $O(n^{3})$ terms, and finally (4.14) has $O(n^{4})$ terms (these orders are true for both: directed and undirected networks). Thus, it is sufficient to show that

[TABLE]

For the proof of (4.18), we note that

[TABLE]

where

[TABLE]

We get

[TABLE]

where, again, (2.10) and continuity of $f_{2}(\theta,t)=-\partial_{\theta^{2}}g(\theta,t)$ has been used. Furthermore, we have that

[TABLE]

by Assumptions (2.11) and (2.18). From (4.19) and (4.23), we obtain (4.18). This shows(4.8).

For the proof of (4.9), we calculate a bound for the expectation of the absolute value of the third derivative of $P_{n}$. With $s=t_{0}+uh$, it holds (recall that $\tau:=\sup_{\theta\in V}\|\theta\|$)

[TABLE]

where (2.6) has been used to get that the order of differentiation and integration can be interchanged and where Fubini could be used because all involved terms are non-negative. The upper bound for the expectation in the integral expression is bounded by Assumptions (2.6) and (2.10). This shows (4.9). ∎

Lemma 4.5**.**

Assume that Fact 4.1, (A3), (A4) (2.7), (2.8), (A6) (2.12), (2.13), (A7) (2.19) hold. It holds that

[TABLE]

With $B_{n}$ from Theorem 2.1, we have

[TABLE]

Proof.

In this proof, we use the shorthand notation K_{h,t_{0}}(s)=\frac{1}{h}K\big{(}\frac{z-t_{0}}{h}\big{)}. We begin with proving (4.25). Denote for vectors $a,b\in\mathbb{R}^{q}$ by $[a,b]$ the connecting line between $a$ and $b$. Note firstly that by a Taylor series application for a random (depending on $X_{n,ij}(s)$) intermediate value $\theta^{*}(s)\in[\theta_{0}(s),\theta_{0,n}]$

[TABLE]

Hence, we obtain

[TABLE]

We decompose (4.27) into two terms by splitting $\theta_{0}(s)-\theta_{0}(t_{0})+\theta_{0}(t_{0})-\theta_{0,n}=(\theta_{0}(s)-\theta_{0}(t_{0}))+(\theta_{0}(t_{0})-\theta_{0,n})$. For the second part we obtain, by using that $\|\theta^{*}(s)\|$ is bounded by $\tau$ (because $\theta_{0,n},\theta_{0}(t_{0})\in V$ by Lemma 4.2 and $\theta_{0}$ is continuous), use also Fubini in the second line and rewrite as a conditional expectation in the last line

[TABLE]

where the last equality holds, because by assumption (2.12) the first factor is $O(h)$, the second factor is uniformly bounded by (2.8) and $\|\theta_{0,n}-\theta_{0}(t_{0})\|=O(h^{2})$ by Proposition 4.2. We now discuss the second term of the split of (4.27). Recall therefore the definitions of $\gamma_{n,ij}(s)$ and $\tau_{n,ij}(\theta,s)$ from Theorem 2.1 and (2.15), respectively. Applying the above and using that $\theta_{0}(s)-\theta_{0}(t_{0})=\theta^{\prime}_{0}(t^{*})(s-t_{0})$ for an appropriate point $t^{*}\in[t_{0},s]$, we obtain

[TABLE]

Hence, we need to prove that (4.31) and (4.32) are $o_{P}(1)$ (these lines individually without the leading $h^{2}$ from the first line) and we are done with the proof. $K$ is supported on $[-1,1]$ and hence $s\in U_{h}:=[t_{0}-h,t_{0}+h]$. Moreover, continuity of $\theta_{0}^{\prime}$ yields $\sup_{s\in U_{h}}\frac{(\theta_{0}(t^{*})-\theta_{0}^{\prime}(t_{0})(t_{0}-s)}{h}\to 0$. Hence, we can show $\eqref{eq:teil1}=o_{P}(1)$ by similar arguments which lead to (4.30). For (4.32) we apply apply Taylor again to get for another intermediate point $\theta^{**}(s)\in[\theta_{0}(s),\theta^{*}(s)]$

[TABLE]

Now arguments are again similar to the ones leading to (4.30), we just have to use the power three part in (2.8) and the fact that $\sup_{s\in U_{h}}\|\theta^{*}(s)-\theta_{0}(s)\|\leq\sup_{s\in U_{h}}\|\theta_{0}(s)-\theta_{0,n}\|$ which converges to zero by continuity of $\theta$ and Proposition 4.2. This concludes the proof of (4.25).

To prove (4.24), we have to show that

[TABLE]

where $r_{n,ij}(s)$ was defined before Assumption (A7). We do this by showing that every component of the left hand side of (4.33) is $o_{P}(1)$, i.e., we replace $r_{n,ij}(s)$ by $r_{n,ij}^{(a)}$ for $a\in\{1,...,q\}$. By an application of Markov’s inequality, this holds if

[TABLE]

We show this similarly as in the proof of Lemma 4.4 by splitting the sum in three sums corresponding to $|\{i,j\}\cap\{k,l\}|=2,1,$ or [math]. The corresponding sums have $O(n^{2}),O(n^{3})$ and $O(n^{4})$ terms, respectively. Before going through these three cases, we note that equations (4.26) and (2.7) imply that $\sup_{u,v\in[-1,1]}\mathbb{E}(r_{n,ij}^{(a)}(t_{0}+uh)r_{n,kl}^{(a)}(t_{0}+vh)|C_{n,ij}(t_{0})=1,C_{n,kl}(t_{0})=1)=O(h^{2})$ for all $a\in\{1,...,q\}$ and for all $(i,j)$ and $(k,l)$. Now we get for the sum over edges with $|\{i,j\}\cap\{k,l\}|=2$ the bound

[TABLE]

For the sum over edges with $|\{i,j\}\cap\{k,l\}|=1$, we get the following bound from (2.13)

[TABLE]

Observing that $\frac{l_{n}h^{3}}{n}=\frac{l_{n}^{3/5}(h^{5})^{3/5}l_{n}^{2/5}}{n}=O(\frac{l_{n}^{2/5}}{n})=O(n^{-1/5}P(C_{n,12}(t_{0})=1)^{2/5})=o(1)$, the bound is of order $o(1).$

By using (2.13) and (2.19), we get the following bound for the sum over edges with $|\{i,j\}\cap\{k,l\}|=0$:

[TABLE]

This concludes the proof of (4.24). ∎

Proposition 4.3**.**

Assume that the assumptions of the Lemmas 4.4 and 4.5 hold. With probability tending to one, the equation $P^{\prime}_{n}(\theta)=0$ (cf. equation (4.2)) has a solution $\widetilde{\theta}_{n}(t_{0})$, which has the property

[TABLE]

To prove this proposition, we will make use of the following theorem, see [5]:

Theorem 4.4**.**

(Newton-Kantorovich Theorem) Let $R(x)=0$ be a system of equations where $R:D_{0}\subseteq\mathbb{R}^{q}\to\mathbb{R}$ is a function defined on $D_{0}$. Let $R$ be differentiable and denote by $R^{\prime}$ its first derivative. Assume that there is an $x_{0}$ such that all expressions in the following statements exist and such that the following statements are true

$||R^{\prime}(x_{0})^{-1}||\leq B$, 2. 2.

$||R^{\prime}(x_{0})^{-1}R(x_{0})||\leq\eta$, 3. 3.

$||R^{\prime}(x)-R^{\prime}(y)||\leq K||x-y||$* for all $x,y\in D_{0}$,* 4. 4.

$r:=BK\eta\leq\frac{1}{2}$* and $\Omega_{*}:=\{x:||x-x_{0}||<2\eta\}\subseteq D_{0}$.*

Then there is $x^{*}\in\Omega_{*}$ with $R(x^{*})=0$ and

[TABLE]

Proof of Proposition 4.3.

We show that $P^{\prime}_{n}(\theta)$ has a root by using Theorem 4.4 with $D_{0}=V$ and $x_{0}=\theta_{0,n}$. Lemma 4.4 gives that $P^{\prime}_{n}(\theta_{0,n})\overset{\mathbb{P}}{\to}0$ and $P^{\prime\prime}_{n}(\theta_{0,n})\overset{\mathbb{P}}{\to}-\Sigma$. Since $\Sigma$ is invertible we also have that the sequence of random variables $B_{n}:=||P_{n}^{\prime\prime}(\theta_{0,n})^{-1}||$ is well-defined (for large $n$) and that it is of order $O_{P}(1)$. Thus we also have $\eta_{n}:=||P^{\prime\prime}_{n}(\theta_{0,n})^{-1}P^{\prime}_{n}(\theta_{0,n})||=o_{P}(1)$. For the Lipschitz continuity of $P_{n}^{\prime\prime}$ we bound the partial derivatives of $P_{n}^{\prime\prime}$ by Lemma 4.4. Hence we conclude that every realization of $P_{n}^{\prime\prime}$ is Lipschitz continuous with (random) Lipschitz constant $K_{n}=O_{P}(1)$. Combining everything, we get that $r_{n}:=B_{n}K_{n}\eta_{n}=o_{P}(1)$. Thus with probability tending to one we have $r_{n}\leq\frac{1}{2}$, and hence the Newton-Kantorovich Theorem tells us that with probability tending to one the equation $P^{\prime}_{n}(\theta)=0$ has a solution $\widetilde{\theta}_{n}(t_{0})\in D_{0}=V$ with the property that

[TABLE]

To prove the asserted rate, we have to investigate $\eta_{n}$ further. We note first that since $P_{n}^{\prime\prime}(\theta_{0,n})^{-1}$ is stochastically bounded, the rate of $\eta_{n}$ is determined by the rate with which $P^{\prime}_{n}(\theta_{0,n})$ converges to zero. To find this rate we observe that every summand of $P_{n}^{\prime}(\theta_{0,n})$ has expectation zero conditionally on $C_{n,ij}(s)=1$:

[TABLE]

by the assumption that $\theta_{0,n}$ maximizes $\theta\mapsto\int_{0}^{T}K\left(\frac{s-t_{0}}{h}\right)g(\theta,s)ds$. So, in $P^{\prime}_{n}(\theta_{0,n}),$ we can subtract $C_{n,ij}(t_{0})\int_{0}^{T}K\left(\frac{s-t_{0}}{h}\right)\partial_{\theta}g(\theta_{0,n},s)ds$ from every summand without changing anything, i.e.,

[TABLE]

By Lemma 4.5, this term is equal to $h^{2}\cdot B_{n}+o_{P}\left(\frac{1}{\sqrt{l_{n}h}}\right)+o_{P}(h^{2})$, which concludes the proof of Proposition 4.3. ∎

Lemma 4.6**.**

Assume that the assumptions of Lemmas 4.4 and 4.5 hold. For $k,l\in\{1,...,q\},$ we have that

[TABLE]

and

[TABLE]

Moreover, it holds that

[TABLE]

Proof.

The proof of (4.34) follows by using similar arguments as in the proof of Lemma 4.4, with $\theta_{0,n}$ replaced by $\theta_{0}(s)$, and with $K$ replaced by $K^{2}$.

For the proof of claim (4.35), we calculate the expectation of the left hand side of (4.35). Because the integrand is positive, we can apply Fubini, and we get that the expectation is equal to

[TABLE]

Here we use (2.10), $\max_{-1\leq u\leq 1}K(u)<\infty$ and (2.9). This shows (4.35).

To see (4.36), we show that

[TABLE]

This then implies (4.36) because of (4.8).

By using exactly the same arguments as in the proof of Lemma 4.5, we obtain

[TABLE]

This gives

[TABLE]

The expectation of the first factor is bounded because of assumptions (2.10) and (2.9). Furthermore, the second term is of order $o_{P}(1)$ by Proposition 4.3. Thus, the product is of order $o_{P}(1)$. This shows (4.37) and concludes the proof of (4.36). ∎

Proposition 4.5**.**

Assume the assumptions of Proposition 4.3 and Lemma 4.6 hold. With probability tending to one, $\partial_{\theta}\ell_{T}(\theta,t_{0})=0$ has a solution $\hat{\theta}_{n}(t_{0})$, and

[TABLE]

Proof of Proposition 4.5.

The proof is based on modifications of arguments used in the asymptotic analysis of parametric counting process models, see e.g. the proof of Theorem VI.1.1 on p.  422 in [2]. Define

[TABLE]

and let $U_{t}^{l}(\theta)$ be defined as $U^{l}(\theta),$ but with $t$ being the upper limit of the integral in (2.2), (i.e., $U^{l}(\theta)=U_{T}^{l}(\theta)$). Furthermore, we write $U(\theta)=(U^{1}(\theta),...,U^{q}(\theta)),$ and the vector $U_{t}(\theta)$ is defined analogously. In the first step of the proof, we will show that

[TABLE]

For the local, square integrable martingale $M_{n,ij}$ defined in (4.1), it holds that $M_{n,ij}$ and $M_{n,i^{\prime}j^{\prime}}$ are orthogonal, meaning that $<M_{n,ij},M_{{n,i^{\prime}j^{\prime}}}>_{t}=0$ if $(i,j)\neq(i^{\prime},j^{\prime})$, i.e. the predictable covariation process is equal to zero. For the predictable variation process of $M_{n,ij},$ we have

[TABLE]

By definition of $\widetilde{\theta}_{n}(t_{0}),$ see the statement of Proposition 4.3, we have that (write $K_{h,t_{0}}(s):=K\left(\frac{s-t_{0}}{h}\right)$)

[TABLE]

So $\widetilde{\theta}_{n}(t_{0})$ was chosen such that the non-martingale part of $\partial_{\theta}\ell(\widetilde{\theta}_{n}(t_{0}),t_{0})$ vanishes. Now, we want to apply Rebolledo’s Martingale Convergence Theorem, see e.g. Theorem II.5.1 in [2]. This theorem implies (4.38), provided a Lindeberg condition (4.35) holds, and

[TABLE]

To verify (4.41), first note that (4.39) and (4.34) imply finiteness of

[TABLE]

with probability tending to one. Note that Lemma 4.6 is formulated with $t=T$, but the integral is finite also for $t<T$ simply because the integrand is non-negative. From now on we assume the above integral is finite. The process

[TABLE]

is a local square integrable martingale, see e.g. Theorem II.3.1 on p.71 in [2]. Since the martingales $M_{n,ij}$ are orthogonal, and by using Lemma 4.6, the predictable covariation satisfies

[TABLE]

This shows (4.41), and concludes the proof of (4.38).

We now show that

[TABLE]

where

[TABLE]

We want to apply the Newton-Kantorovich Theorem 4.4 with $R(\theta):=R_{n}(\theta):=\frac{1}{l_{n}h}U(\theta)$, $D_{0}=V$ and $x_{0}:=\widetilde{\theta}_{n}(t_{0})$. To this end, define

[TABLE]

From Lemma 4.6, we know that $P^{\prime\prime}_{n}(\widetilde{\theta}_{n}(t_{0}))$ converges and is invertible for $n$ large enough, and thus $B_{n}=O_{P}(1)$. Now let

[TABLE]

Results (4.36) and (4.38) imply that $\eta_{n}=o_{P}(1)$. Next, notice that $P^{\prime\prime}_{n}$ has a Lipschitz constant $K_{n}$ that is bounded by the maximum of the third derivative of $P_{n}$. According to (4.9), this maximum is bounded, and we obtain $K_{n}=O_{P}(1)$. Hence, $r_{n}=B_{n}K_{n}\eta_{n}=o_{P}(1)$. Now, Theorem 4.4 implies that, with probability converging to one, there is $\hat{\theta}_{n}(t_{0})$ such that $U(\hat{\theta}_{n}(t_{0}))=0$ and

[TABLE]

To obtain the asymptotic distribution of $\hat{\theta}_{n}(t_{0}),$ we note that, by (4.38) and (4.36),

[TABLE]

Thus it holds that $\sqrt{l_{n}h}\cdot Z_{n}=O_{P}(1),$ and as a consequence we get $\sqrt{l_{n}h}\cdot\eta_{n}=O_{P}(1)$. Using the second statement of the Newton-Kantorovich Theorem 4.4, we obtain

[TABLE]

Thus $\sqrt{l_{n}h}\cdot(\widetilde{\theta}_{n}(t_{0})-\hat{\theta}_{n}(t_{0}))$ and $\sqrt{l_{n}h}\cdot Z_{n}$ have the same limit distribution. Because of (4.43) this implies the statement of the proposition. ∎

Proof of Theorem 2.1 Combining Propositions 4.3 and 4.5, and applying Slutzky’s Lemma, we obtain by the assumptions on the bandwidth $h$ in (A2)

[TABLE]

With Proposition 4.2 this gives (2.20).

5 Acknowledgement

The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant INST 35/1134-1 FUGG. Research of Alexander Kreiß and Enno Mammen was supported by Deutsche Forschungsgemeinschaft through the Research Training Group RTG 1953. Research of W. Polonik has been supported by the National Science Foundation under Grant No. DMS 1713108.

We would like to thank Michael Gertz and Andreas Spitz for many helpful discussions and suggestions. In particular, we are very grateful to Andreas Spitz for spending a lot of time in helping us processing and preparing the data.

6 Appendix

6.1 Simulations of degree distributions, cluster coefficients and diameters.

Here we report additional simulations of degree distributions, cluster coefficients and diameters. In Section 3, we have presented results for the degree distribution of networks based on the Washington DC bikeshare activity on 7th December 2012. In this section, we will consider the days 18th April 2014 and 10th July 2015, and also compare diameters and clustering coefficients of the simulated and observed networks. As above, using the corresponding estimated parameter value for each of these days, we compute 3840 predictions and compared them with the observed values. The diameter of a network is the longest among the shortest path between two vertices in the network. Typically, in observed networks the diameter is much smaller than the number of vertices (cf. [15]). The clustering coefficient is the number of complete triangles (triples of vertices which are completely connected) divided by the number of incomplete triangles (triples of vertices with at least two edges). Note that every complete triangle is also incomplete, hence the clustering coefficient is between zero and one. The clustering coefficient can be understood as the empirical probability that vertices are connected given that there is a third vertex to which both are connected. It has been reported (cf. [15]), that in observed networks this number is usually significantly higher than in an Erdös-Rényi network, where the presence of edges are i.i.d. random variables.

Our question here is, Does a network which was simulated by our model look like the observed network? or in other words Could one believe that the observed network is a realization of our model?. To answer this, we consider the three network characteristics mentioned above, and empirically and visually compare the simulated results to the observed data. The heuristic justification underlying this approach is, that, if considered jointly, these three characteristics are able to discriminate between a range of different types of networks (see also [15, 41])

We start by presenting results for diameter and clustering coefficient on 7th of December 2012. As described in Section 3, where the degree distribution was discussed, we divide the edges between bike stations in six regimes by considering tour frequencies between the stations on the day. Figure 14 shows the histograms of the simulated diameter in the different regimes. We see that, in 14(e) (as before in Figure 5(e)), the simulation and the reality appear to coincide nicely. In other words, for a moderate number of tours our model seems to fit well. It is interesting to note that our model performs differently in the different regimes suggesting that edges with different activity have to be modeled differently. Finally, in Figure 15, we see the histograms of the simulated clustering coefficients. The true value in the corresponding regime is shown in the titles of the plots. Overall, the performance appears reasonable. In particular, in Figure 15(d) the histogram is nicely centered around the true value. Interestingly, the performance in the fifth regime ($l_{1}=5$ and $l_{2}=12$), shown in Figure 15(e), is not as good as the others. One explanation for this might be that here different covariates are needed.

In Figure 22(a) we see one simulated graph compared to the true graph. The color of the edges determine how many tours happened relative the the other edges: The lowest 25% of the edges are colored green, the next 25% yellow, then orange and the highest 25% of edges are colored red. Due to the integral value of the activity it is not the case that exactly 25% of the edges are green and so on. The size of the vertices is relative to their degree. We see that the model is able to find the important (i.e. high degree) vertices. For the edges we see that some red edges are at wrong places. But generally the vertices with high profile edges are recognized. The remaining graphs in Figure 22 show the same comparison for the two other dates under consideration. And we see that the results are similar.

Figures 16 till 21 show the results of the corresponding simulations for the other two dates. Overall the results are similar. It should be pointed out that even though the model is not able to reproduce every feature perfectly accurate, the simulated networks are still very close to the true observation. This becomes more obvious if we remind ourselves that only six parameters were used.

6.2 Bandwidth choice

Under our assumptions that the covariates stay constant over the day, it makes sense to consider only integral bandwidth lengths (for us one day has length one). In order to choose the bandwidth, we apply a one-sided cross validation (cf. [11, 27]) approach which was shortly motivated in Section 3 and which we now describe in detail.

Let $K$ and $L$ be kernels fulfilling the assumptions in the paper and denote by $\hat{\theta}_{K}(t_{0})$ and $\hat{\theta}_{L}(t_{0})$ the maximum likelihood estimators using $K$ and $L$ respectively. Then, by Theorem 2.1, we get that asymptotically the bias and the variance of the estimators can be written as

[TABLE]

where the constants $C_{1}$ and $C_{2}$ depend on the true parameter curve $\theta_{0}$ and the time $t_{0}$ but not on the kernel. Hence, the corresponding expressions for $\hat{\theta}_{L}(t_{0})$ can be found, just by replacing every $K$ with an $L$. The decomposition of the asymptotic mean squared error in squared bias plus variance yields the following asymptotically optimal bandwidths $h_{K}$ and $h_{L}$ , minimizing the asymptotic mean squared error:

[TABLE]

Again, the corresponding expression for $h_{L}$ can be found by replacing every $K$ by $L$. So the following formula, known from kernel estimation, holds also true in our setting

[TABLE]

This means that knowledge of the bandwidth minimizing the mean squared error for kernel $L,$ implies knowledge of the bandwidth minimizing the mean squared error using kernel $K$. Ultimately, we use a triangular kernel $K(u)=(1+u)\mathbbm{1}_{[-1,0)}(u)+(1-u)\mathbbm{1}_{[0,1]}(u)$. In order to find the bandwidth $h_{K}$ for this kernel, we want to apply cross-validation. As proposed in [11] one-sided cross validation is an attractive method for the case of time series data. One-sided here means that we apply cross validation to a kernel $L$ which is only supported on the past $[-1,0]$. In order to avoid a bias, we use the one-sided kernel together with local linear approximation. This following heuristic derivations motivates this choice.

Firstly, in our regular maximum likelihood setting, we maximize, over $\mu\in\Theta,$ the expression

[TABLE]

Deriving this expression with respect to $\mu$, setting the derivative equal to zero, and rearranging terms, yields (to save space we use here a fraction, although the denominator is a matrix)

[TABLE]

Using the notation $y_{1}:=\mathbb{E}(X_{n,ij}(t_{0})e^{\theta_{0}(t_{0})^{T}X_{n,ij}(t_{0})}|C_{n,ij}(t_{0})=1)\cdot\mathbb{P}(C_{n,ij}(t_{0})=1)$ and $y_{2}:=\mathbb{E}(X_{n,ij}(t_{0})X_{n,ij}(t_{0})^{T}e^{\theta_{0}(t_{0})^{T}X_{n,ij}(t_{0})}|C_{n,ij}(t_{0})=1)\cdot\mathbb{P}(C_{n,ij}(t_{0})=1),$ we obtain the approximation

[TABLE]

Now define the local linear estimator $\hat{\theta}_{LC,K}(t_{0}),$ with respect to a kernel $K,$ as the value of $\mu_{0}$ maximizing the following expression over $(\mu_{0},\mu_{1})\in\Theta^{2}:$

[TABLE]

Using the same approximations as in the usual kernel estimation setting, and deriving the resulting approximate likelihood, we obtain

[TABLE]

where $M_{k}:=\int_{-1}^{1}u^{k}K(u)du$. The previous computations were just a heuristic. But nevertheless, the similarity between the previous display and (6.2) suggests that the local linear estimator $\hat{\theta}_{LC,K}$ using the kernel $K$ is actually just a regular kernel estimator $\hat{\theta}_{L}$ with kernel function

[TABLE]

This aligns with results about kernel estimation, as, for example, stated in [27]. It can be easily computed that this new kernel is of order one, i.e., $\int uL(u)du=0$, even though the original kernel was not. Hence, knowledge of the optimal bandwidth for the local linear estimator using the kernel $K$ implies knowledge of the optimal bandwidth for any other order one kernel by means of (6.1). Taking the same route as in [11], the selector for the bandwidth $\hat{h}_{K}$ for the triangular kernel $K$ is the following: Let $K^{*}(u):=2K(u)\mathbbm{1}_{[-1,0]}(u)$ denote the one sided version of $K.$

Find a bandwidth $\hat{h}_{L}$ for the local linear estimator $\hat{\theta}_{LC,K^{*}}$ based on the kernel $K^{*}$ via cross validation (since we use a one-sided kernel, this step is also called one-sided cross-validation. We will make it more precise later). 2. 2.

Compute $\hat{h}$ by using (6.1) with $L$ defined as in (6.3) but with $K$ replaced by $K^{*}$.

For the one-sided cross-validation in step 1, we minimize, in our bike share data analysis, the following function in $h$

[TABLE]

where $k_{T}$ was the number of weeks (recall that we assume the covariates to remain constant over a day, and that we only consider Fridays), i.e., $k$ refers to the $k$-th Friday in the dataset. $L(k)$ refers to the set of pairs $(i,j)$ between which there was a bike tour on Friday $k$, $X_{i,j}(k)$ is the true number of bike tours observed between $i$ and $j$ on Friday $k$. Finally, $\hat{\theta}_{LC,K^{*}}^{(-k)}(k)$ is the local linear estimator with respect to the kernel $K^{*}$ based on all but the $k$-th Friday. Since $K^{*}$ is left-sided, this really means the estimator is based on Fridays $0,...,k-1,$ and hence the term one-sided cross-validation. The intensities are the theoretical values of the expectation of the number of bike tours if the model is correct. So we compute the squared difference with the true number of bike rides and divide by the estimated intensity, where we only take the non-censored edges into account.

In Section 3, we had displayed results for different bandwidths $h$ of (6.4) in Figure 6. The prediction error of the fit decreases, until the bandwidth is equal to 23. Afterwards the prediction error stays roughly the same and starts to increase when the bandwidth reaches a full year (52 weeks). This may be explained by a periodicity with a period of approximately one year. If one uses 23 as minimal value we get as asymptotic optimal bandwidth 12 which is approximately 23 divided by $\rho$. Here, following Step 2 of the above described procedure, we use that $\rho$ is approximately equal to 1.82 for triangular kernels.