Estimating Network Structure from Incomplete Event Data
Abstract
Multivariate Bernoulli autoregressive (BAR) processes model time series of events in which the likelihood of current events is determined by the times and locations of past events. These processes can be used to model nonlinear dynamical systems corresponding to criminal activity, responses of patients to different medical treatment plans, opinion dynamics across social networks, epidemic spread, and more. Past work examines this problem under the assumption that the event data is complete, but in many cases only a fraction of events are observed. Incomplete observations pose a significant challenge in this setting because the unobserved events still govern the underlying dynamical system. In this work, we develop a novel approach to estimating the parameters of a BAR process in the presence of unobserved events via an unbiased estimator of the complete data loglikelihood function. We propose a computationally efficient estimation algorithm which approximates this estimator via Taylor series truncation and establish theoretical results for both the statistical error and optimization error of our algorithm. We further justify our approach by testing our method on both simulated data and a real data set consisting of crimes recorded by the city of Chicago. ^{†}^{†}This research is funded by the National GeospatialIntelligence Agency’s Academic Research Program (Award No. HM04761712003, Project Title: Inference for Autoregressive Point Processes). Approved for public release, 18935. This work was also supported by NSF CCF1418976, NIH 1 U54 AI11792401, NSF CCF0353079, ARO W911NF1710357, and AFOSR FA95501810166.
1 Introduction
Discrete event data arises in a variety of forms including crimes, health events, neural firings, and social media posts. Frequently each event can be associated with a node in a network, and practitioners aim to learn the relationships between the nodes in the network from the event data. For example, one might observe a sequence of crimes associated with different gangs and seek to learn which crimes are most likely to spark retaliatory violence from rival gangs.
Such problems have attracted widespread interest in recent decades and a variety of point process models have been proposed to model such data. The Hawkes process ([13]) is a popular continuous time framework which has been applied in a number of different contexts (e.g., [34, 33, 8, 31]). In addition, other works works have used a discrete time framework to model time series event data (e.g., [17, 10, 3]).
A central assumption of many of these works is that all the events are observed. However, in many cases we may observe only a subset of the events at random. For example, while point process models have been widely used to model crime incidence ([23, 24, 18]), frequently one only has access to reported crime data. For many crimes the true number of incidents can be substantially higher. The gap between the reported and true crime rates is referred to as “the dark figure of crime” by researchers in Sociology and Criminology who have studied this issue extensively ([6, 15]). Unobserved events also pose a challenge in inference from Electronic Health Record (EHR) data which can be incomplete for a number of different reasons ([30, 29]).
The unobserved events still play a role in the dynamical system governing the time series, making network estimation with incomplete data particularly challenging. In this paper, we contribute to the growing literature on modeling in the presence of unobserved events by proposing a novel method for estimating networks when we only observe a subset of the true events.
1.1 Problem Formulation
Many point process models of time series of discrete events are temporally discretized either because event times were discretized during the data collection process or for computational reasons. In such contexts, the temporal discretization is typically such that either one or zero events are observed in each discrete time block for each network node. With this in mind, we model the true but unobserved observations using a Bernoulli autoregressive process:
(1.1) 
Here is a vector indicating whether events occurred in each of the nodes during time period . The vector is a constant bias term, and the matrix is the weighted adjacency matrix associated with the network we wish to estimate. We assume that each row of lies in the ball of radius , which we denote . We generally consider a case where is sparse and the magnitude of all its entries are bounded, so that is a universal constant which is independent of .
We observe , a corrupted version of (1.1) where only a fraction of events are observed as follows:
(1.2) 
Here denotes the Hadamard product and is a vector where each entry is independently drawn to be one with probability and zero with probability .
Our analysis of (1.1) and (1.2) can be naturally extended to several more complex variants. Instead of assuming each is observed with probability , we can assume events from each node are observed with a unique probability . We consider only a first order AR(1) process but our framework can be extended to incorporate more sophisticated types of memory as in [22]. We omit discussion of these extensions in the interest of clarity.
2 Related Work
Corrupted or missing data in highdimensional data sets is a problem which appears in a number of different domains and has attracted widespread interest over the past few decades (see [11] for an applicationfocused overview). Our focus is on a particular type of corrupted data: partially observed sequences of discrete events. In recent years researchers have started to focus on this problem ([32, 28, 16]). The prior works of which we are aware use a Hawkes process framework and assume knowledge of the time periods when the data is corrupted. In the context of (1.2) this amounts to knowledge of . Our method can operate in a setting where the researcher cannot differentiate between periods when no event actually occurred, and when events potentially occurred but were not recorded. Moreover, because we use a discretetime framework, we are able to derive sample complexity bounds for the estimation procedure proposed in Section 3. Our theoretical results complement the empirical nature of much of the past work in this area.
This paper is also related a variety of works on regularized estimation in highdimensional statistics, including [5, 27, 4] and [14]. Many of these works have derived sample complexity guarantees using linear models, and some of these results have been extended to autoregressive generalized linear models ([25, 12]). Another line of research ([20, 1, 21, 19, 25]) has formalized a notion of Restricted Strong Convexity (RSC) which we leverage in Section 4.2. While many loss functions of interest in highdimensional statistics are not strongly convex, these works have shown that they frequently satisfy an RSC condition which is sufficient to provide statistical and optimization error guarantees. The main technical challenges in our setting lie in establishing results similar to these RSC conditions.
2.1 Missing Data in a HighDimensional Linear Model
[20] straddles the missing data literature and highdimensional statistics literature. The authors consider a missing data linear model
(2.1) 
where and one observes pairs and aims to estimate . The authors propose minimizing a loss function of the observed data which satisfies the property
for any . Here
denotes the classical Lasso loss function with the unobserved data and regularization parameter . In other words, the missing data loss function is an unbiased estimator for the full data Lasso loss function we would ideally like to minimize. This idea motivates our construction of a loss function for the observed process (1.2).
Our problem can be viewed as an extension of [20] to autoregressive GLM models without knowledge of .^{1}^{1}1Note that [20] does consider AR processes, but in a different context from our setting. Specifically, we wish to estimate the AR process parameters, where as they consider a special case of (2.1) where but where is known and one aims to estimate . In particular, we cannot distinguish events that were missed ( and ) from correctly observed periods with no events (). [20] are able prove sample complexity bounds as well as optimization bounds which are consistent with the highdimensional statistical literature in that they scale with rather than the dimension of . We are able to prove analogous bounds for our estimator in Section 4.
2.2 Contributions
This paper makes the following contributions.

We propose a novel method for estimating a network when only a subset of the true events are observed. In contrast to previous work, our methods do not rely on knowledge of when the data is potentially missing. Our procedure uses Taylor series approximations to an unbiased loss function, and we show that these approximations have controlled bias and lead to accurate and efficient estimators.

We prove bounds on both the statistical error and optimization error of our proposed estimation method. The results hinge on showing that our loss function satisfies a restricted strong convexity (RSC) condition. Past work on linear inverse problems with corrupted designs also establish RSC conditions, but these conditions do not carry over to the autoregressive GLM setting.

We demonstrate the effectiveness of our methodology on both simulated data and real crime data.
3 Proposed Estimation Procedure
Given the full data , the negative loglikelihood function is decomposable in the rows of . In other words, if
then where is the th row of and denote the loss function restricted to a specific row. Throughout the paper we slightly abuse notation and let refer to the entire loss function when is a matrix, and let refer to the loss function for a specific row when is a row vector. The loss function for the th row takes the form
where is the partition function for the Bernoulli GLM.
We do not have access to and instead we aim to estimate using the corrupted data . As discussed in Section 2.1, our strategy will be to construct a loss function of which is an unbiased estimator for . In other words, we want to find some function such that for any ,
(3.1) 
In contrast to the Gaussian case discussed in Section 2.1, the Bernoulli partition function is not a polynomial and constructing a function satisfying (3.1) directly is challenging. We adopt a strategy of computing unbiased approximations to truncated Taylor series expansions of and arriving at as a limit of such approximations.
To do this, we first rewrite using its Taylor series expansion around zero
The constant factor does not effect estimation in any way so we ignore it for the remainder of our discussion in the interest of simplicity. We let denote the degree Taylor truncation to . The are binary vectors and we assume each row is sparse, so will not be too far from zero. Thus it is reasonable to hope that for small , is a good approximation for whenever . We bound the approximation error in Lemma B.3 in the supplement.
We now consider the problem of constructing a function such that
(3.2) 
Once we construct we can estimate the th row of by attempting to solve
Key question we need to address with this approach include (a) can we (approximately) solve this optimization problem efficiently? (b) will the solution to this optimization problem be robust to initialization? (c) will it be a strong estimate of the ground truth?
3.1 Definition of
We first derive an unbiased estimator of the degreetwo Taylor series expansion .
Note that there are straightforward unbiased estimates of the first two terms:
(3.3) 
For the third term, , note that
3.2 HigherOrder Expansions
The construction of in the previous section suggests a general strategy for constructing satisfying (3.2). Take any monomial
depending on the counts in nodes during time period . Wherever this monomial appears in , our unbiased loss function will have a term
scaled by where denotes the number of unique terms in the monomial. For example, in Equation (3.3) each degree two monomial was unique so we scaled everything by . However, in (3.5) some of the degree two monomials had repeated terms and so they were scaled by . In order to formalize these ideas and generalize our estimator to , we first need to introduce additional notation.
3.2.1 Notation
Let denote the set of all monomials of degree . We represent an element as a list containing elements. An element in the list corresponds to the index of a term in the monomial. For an example, the monomial can be represented as the list .
For a polynomial function we let denote the coefficient of the monomial in . Finally we define the order of a list to denote the number of unique elements in the list, so whereas .
Example
Consider the function . We can decompose as
where with corresponding coefficients , and .
Using this notation we can write
3.2.2 Definition of
The degree likelihood is constructed as follows
Recall that denotes the number of unique terms in the monomial . In other words, in we adjust by scaling each monomial according the the number of unique terms rather than the number of overall terms. This definition clearly satisfies (3.2). We show in the supplement that if and then converges uniformly on to a function we denote as . Extending this loss function on individual rows to an entire matrix, we can define . An additional technical discussion in the supplement guarantees that actually satisfies the desired property in Equation (3.1).
3.3 Proposed Optimization
In practice we can only compute for finite . To estimate we consider the following constrained optimization:
(3.6) 
where is an estimate of the missingness parameter and
In general, is not a convex function. However, we show in Section 4.2 that under certain assumptions all stationary points of (3.6) must lie near . Thus we can approximately solve (3.6) via a simple projected gradient descent algorithm.
In order to apply our algorithm it is necessary to have an estimate of the frequency of missed data. In many cases one may have prior knowledge available. For example, social scientists have attempted to quantify the frequency of unreported crimes ([15, 26]). Moreover, a simulation study in Section 5 suggests our strategy is robust to misspecification of .
4 Learning Rates
In this section we answer several questions which naturally arise with our proposed estimation procedure in Equation (3.6). In particular, we’d like to know:

If we can find a solution to Equation (3.6), will it be a good approximation to ?

The loss function is not convex. Is it possible to actually find a minimizer to (3.6), or an approximation to it?
In Theorem 4.1 we answer the first question affirmatively. In Theorem 4.2 we show that all stationary points of (3.6) are near one another, so that simple gradientbased methods will give us a good solution.
Throughout this section we assume and . All results in the section apply for the loss functions for . In the case we recover the idealized loss function .
We use to mean and to mean where is a universal constant. Define and .
4.1 Statistical Error
The following theorem controls the statistical error of our proposed estimator.
Theorem 4.1 (Accuracy of ).
Suppose
where . Then
for with probability at least .
The two terms in the upper bound of Theorem 4.2 have a natural interpretation. The first represents the error for the idealized estimator , while the second represents the error due to the Taylor series truncation. Our error scales as which is reasonable in the context of our algorithm because is only welldefined when (see Remark 2 in the supplement). An interesting open question which arises from Theorem 4.1 is whether the process described in (1.1) and (1.2) is unidentifiable for or something specific to our methodology fails for below this threshold.
The proof of Theorem 4.1 uses ideas from the analysis of highdimensional GLMs ([12, 22]) as well as ideas from the analysis of missing data in the linear model [20] and Gaussian linear autoregressive processes [14]. The key technical challenge in the proof lies in controlling the gradient of the error term . This is done in Lemmas B.2B.6 in the supplement.
4.2 Optimization Error
We next focus on the optimization aspects of Equation (3.6). Our loss function is nonconvex, so at first glance it may appear to be a poor proxy loss function to optimize. However, a body of research (see [1, 21, 20, 19]) has studied loss functions satisfying a properties known as restricted strong convexity (RSC) and restricted smoothness (RSM). These works have shown that under certain conditions, the optimization of nonconvex loss functions may be tractable. The formal definitions of the RSC and RSM conditions we use are as follows.
Definition 1 (Restricted Strong Convexity).
Let denote the first order Taylor approximation to a loss function centered at . A loss function satisfies the RSC condition with parameters if
for all .
Definition 2 (Restricted Smoothness).
A loss function satisfies the RSM condition with parameters if
for all .
We are able to show these conditions are satisfied for and . This in turn gives the following result. As in Theorem 4.1 we assume .
Theorem 4.2.
Suppose and for at least rows of where is a universal constant. Let be any stationary point of where . Then
with probability at least for and .
As in Theorem 4.1 the first term in our bound can be interpreted as the error for the idealized estimator while the second term can be thought of as the error due to the Taylor series truncation. The assumption that for at least rows of says that at least a constant fraction of nodes are influenced by other nodes in the network. This assumption allows us to state Theorem 4.2 in terms of  the support of . In extreme cases where almost all nodes in the network fire independently of the other nodes it is possible for the optimization error to have a slower scaling than .
The RSC and RSM conditions are closely related to ideas used in our statistical error bounds in Theorem 4.1. Lemma D.2 shows that the conditions are satisfied for on the order of which leads to an overall optimization error bound of the same order. This is a slower convergence rate than in the linear case; whether stronger rates can be obtained in the autoregressive GLM setting is an open question.
The proof of Theorem 4.2 proceeds as follows. We first establish that the RSC/RSM conditions hold for reasonable constants in Lemma D.2. This proofs relies on the technical machinery built up in Lemmas B.2B.6. We can then combine our RSC/RSM results with Theorem 2 in [1] to conclude that all stationary points of lie in a small neighborhood of with high probability.
5 Simulations
In this section we evaluate the proposed method on synthetic data. We generate with nonzero entries with locations chosen at random. Each nonzero entry is chosen uniformly in and . We then generate a “true" data set and an “observed” data set according (1.1) and (1.2) with . We perform projected gradient descent with a random initialization and show a median of 50 trials.
Figure LABEL:figure:figure2 shows mean squared error (MSE) vs for (top) and (bottom). Our method is shown in red. It uses the loss function on the partially observed data . Our method is compared to the loss function using both the full data (i.e., an oracle estimator with access to the missing data) and the partially observed data (i.e. a naive estimator that ignores the possibility of missing data). As expected, with access to the complete data one can get a more accurate estimate of than either method using the partially observed data. However, our method outperforms the full data likelihood when given the partially observed data. In particular, note that the accuracy for the full data likelihood stalls after some time due to the inherent bias in using the corrupted data on the true data likelihood. In contrast our unbiased method continues to converge to the solution, as suggested by the results in Section 4. Finally, observe that for large there is little variation between trials when using even though each trial was initialized randomly. This agrees Theorem 4.2 which states that all stationary points of lie near one another.
In practical applications one may have strong reason to believe some events are unobserved, but pinning down a precise estimate of the missingness parameter might be unrealistic. Therefore it is important to see how our algorithm performs as a function of the misspecification . We examine this in Figure LABEL:fig:robustness. We generate data as in the previous section but with . We then apply our algorithm with the loss function and varying values of .
Figure LABEL:fig:robustness shows that our method is highly robust to small misspecification of the missingness parameter . Interestingly, underestimating by more than 10% leads to poor results but our method is still robust to overestimation of over 10%. This suggests there is value in applying our techniques with a conservative estimate of the amount of missed data, even when one has only a rough estimate of the frequency of missed events.
As final experiment we measure how MSE varies as a function of the Taylor series truncation level . Calculating takes a significant amount of time for highdimensional problems, so we randomly generate with nonzero entries compared to in the previous simulations and run trials. We set .
In Figure LABEL:figure:figure1 we show MSE as a function of for the full data loss function at three different truncation levels: and . Recall that and has no odd degree terms other than , so and . We see that the second and fourth degree truncations perform essentially the same as the full data likelihood. We also plot MSE as a function of for the truncated missing data loss functions and . As expected, using the full data gives stronger results than the partially observed data. We again see that the second and fourth order truncations perform nearly the same. The sample standard deviations are also similar  e.g. when the standard deviations of and are and respectively while the standard deviations for and are and . The similarity between the second and fourth order truncation levels suggests that choosing one of these truncation levels will give us a strong approximation to . Since takes significantly longer to compute, we use the second order approximation in the first two experiments.
6 Chicago Crime Data
This section studies a data set consisting of crimes committed in Chicago between January 2001 and July 2018 ([9]). Point process models have been applied to this data set in the past ([18]). In a missing data setting, in order to validate our model it is important to have a “ground truth” data set. For this reason we limit our study to homicides within the data set. For other crimes recorded data is known to be incomplete ([15, 6]), but we assume that nearly every murder is observed. This allows us to create an incomplete data set by removing murders randomly while still maintaining a ground truth data set for validation.
The city is divided into 77 community areas and the community area where each murder occurred is recorded. The majority of these community areas experience few murders so we focus on the nine areas with the most murders since 2001. These areas form two clusters: one on the west side of the city and another on the south side. We discretize the events using one week bins, so if a murder occurred in community area during week and otherwise. This gives a data matrix which we divide into a train set containing the first 600 weeks in the period, and a test set containing the final 318 weeks. We then create an incomplete data set where contains independent realizations of a Bernoulli random variable with mean .
We learn parameters and using the training set and the full data likelihood . We also learn parameters and using the incomplete train data and the missing data likelihood for various values of .
We compare the loglikelihood of these parameters on the test set . The results are shown in Figure 1. The missing data estimates perform nearly as well as the full data estimate when is close to the true value of . Note that closely approximates the full data likelihood and the hold out likelihood is substantially worse for compared to for close to . In other words, ignoring the missing data entirely gives a weaker estimate than applying the techniques this paper introduces, even when is not a perfect estimate of . Finally, observe that is more robust to misspecification when compared to when . This is a trend which also appears in Figure LABEL:fig:robustness and suggests there is value in using conservative estimates of the amount of missed data in practical applications.
Given estimates of and we can use density propagation to predict the likelihood of homicides during week based on observed homicides up to week . We do this for learned from the incomplete data with as well as learned from but with , which corresponds to assuming there is no missing data. We use particle filtering to construct estimates
and
These probabilities correspond to the likelihood of homicides during the th week based on the observations over the first weeks. We construct such estimates for each week in the week test set. As expected assigns higher likelihoods of homicides, with expected homicides in total compared to for . As a naive method of correcting for missing data, we divide by a constant scaling factor of and report these likelihoods below; by doing this, we ensure that both predictions yield similar average numbers of homicides, so differences in performance between the proposed and naive estimator are not due to a simple difference in averages, but rather because the proposed method is capturing the underlying dynamics of the system.
Figure 2 displays these likelihoods for Community Area 25 (Austin) which has the largest number of homicides recorded during the test period. We use Gaussian smoothing to help visualize the trends. The top panel shows the predicted probability of events using (in red) and the scaled predicted probability of events using (in blue). The bottom panel shows the actual and incompletely observed events during the test period. The true events generally peak at times in which the predicted events for peak. For example, both the predicted event and true event charts have peaks around weeks and . In contrast, the predicted events for are nearly constant over time. Since it does not account for the missing data (except via a uniform scaling factor), the network is not able to capture the dynamics of the process and so it cannot predict events with as much precision as .
7 Conclusion
We propose a novel estimator for Bernoulli autoregressive models which accounts for partiallyobserved event data. This model can be used in a variety of contexts in which observed discrete event data exhibits autoregressive structure. We provide mean squared error bounds which suggest that our method can accurately capture a network’s structure in the presence of missing events. Simulations and a real data experiment show that our method yields significant improvement compared with ignoring the missed events and that our method is robust to misspecification of the proportion of missed events. The framework described in this paper suggests a strategy for addressing regression problems with corrupted data in other GLMs, although further work is needed to extend our theoretical analysis beyond binary observations.
References
 [1] A Agarwal, S Negahban, and M. Wainwright. Fast global convergence of gradient methods for highdimensional statistical recovery. The Annals of Statistics, 40(5):2452–2482, 2012.
 [2] H. Alzer. Sharp bounds for the bernoulli numbers. Archiv der Mathematik, 74(3):207–211, 2000.
 [3] Y. Bao, C. Kwong, P. Peissig, D. Page, and R. Willett. Point process modeling of adverse drug reactions with longitudinal observational data. In Proc. Machine Learning and Healthcare, 2017.
 [4] S. Basu and G. Michailidis. Regularized estimation in sparse highdimensional time series models. Annals of Statistics, 43(4):1535–1567, 2015.
 [5] P. Bickel, Y. Ritov, and A. Tsybakov. Simultaneous analysis of Lasso and Dantzig selector. Annals of Statistics, 37(4):1705–1732, 2009.
 [6] A. Biderman. Understanding Crime Incidence Statistics. SpringerVerlag, New York, first edition, 1991.
 [7] L. Carlitz. Eulerian numbers and polynomials. Mathematics Magazine, 32(5):247–260, 1959.
 [8] S. Chen, A. Shojaie, E. SheaBrown, and D. Witten. The multivariate Hawkes process in high dimensions: Beyond mutual excitation. arXiv preprint arXiv:1707.04928, 2017.
 [9] City of Chicago. Crimes  2001 to present, 2018. https://data.cityofchicago.org.
 [10] A. Fletcher and S. Rangan. Scalable inference for neuronal connectivity from calcium imaging. In NIPS, 2014. arXiv:1409.0289.
 [11] J. Graham. Missing data analysis: Making it work in the real world. Annual Review of Psychology, 60:549–576, 2009.
 [12] E. C. Hall, G. Raskutti, and R. Willett. Inference of highdimensional autoregressive generalized linear models. arXiv preprint arXiv:1605.02693, 2016.
 [13] A. G. Hawkes. Point spectra of some selfexciting and mutuallyexciting point processes. Journal of the Royal Statistical Society. Series B (Methodological), 58:83–90, 1971.
 [14] A. Jalali and R. Willett. Missing data in sparse transition matrix estimation for subgaussian vector autoregressive processes. In American Control Conference, 2018. arXiv:1802.09511.
 [15] L. Langton, M. Berzofsky, C. Kerbs, and H. SmileyMcDonald. Victimizations not reported to the police, 20062010. Technical report, Bureau of Justice Statistics, 2012.
 [16] T. Le. A multivariate Hawkes process with gaps in observations. IEEE Transactions on Information Theory, 63(3):1800–1811, 2018.
 [17] S. Linderman, R. Adams, and J. Pillow. Bayesian latent structure discovery from multineuron recordings. In Advances in neural information processing systems, 2016.
 [18] S. W. Linderman and R. P. Adams. Discovering latent network structure in point process data. In ICML, 2014. arXiv:1402.0914.
 [19] PoLing Loh. Statistical consistency and asymptotic normality for highdimensional robust Mestimators. Annals of Statistics, 42(5):866–896, 2017.
 [20] PoLing Loh and Martin J Wainwright. Highdimensional regression with noisy and missing data: Provable guarantees with nonconvexity. The Annals of Statistics, 40:1637–1664, 2012.
 [21] PoLing Loh and Martin J Wainwright. Regularized mestimators with nonconvexity: Statistical and algorithmic theory for local optima. Journal of Machine Learning Research, 16:559–616, 2015.
 [22] B. Mark, G. Raskutti, and R. Willett. Network estimation from point process data. arXiv preprint arXiv:1802.09511, 2018.
 [23] G. Mohler. Marked point process hotspot maps for homicide and gun crime prediction in chicago. International Journal of Forecasting, 30(3):491–497, 2014.
 [24] G Mohler, M Short, Sean Malinowski, Mark Johnson, G Tita, Andrea Bertozzi, and P. Brantingham. Randomized controlled field trials of predictive policing. Journal of the American Statistical Association, 110(512):1399–1411, 2016.
 [25] Sahand Negahban, Pradeep Ravikumar, Martin Wainwright, and Bin Yu. A unified framework for highdimensional analysis of Mestimators with decomposable regularizers. Statistical Science, 27(4):538–557, 2010.
 [26] T. Palermo, J. Bleck, and A. Peterman. Tip of the iceberg: Reporting and genderbased violence in developing countries. American Journal of Epidemiology, 179(5):602–612, 2014.
 [27] G. Raskutti, M. J. Wainwright, and B. Yu. Restricted eigenvalue conditions for correlated Gaussian designs. Journal of Machine Learning Research, 11:2241–2259, 2010.
 [28] C. Shelton, Z. Qin, and C. Shetty. Hawkes process inference with missing data. In Proceedings of the ThirtySecond AAAI Conference on Artificial Intelligence, 2018.
 [29] NG. Weiskopf and Weng C. Methods and dimensions of electronic health record data quality assessment: Enabling reuse for clinical research. Journal of the American Medical Informatics Association, 20(1):144–151, 2013.
 [30] B. Wells, K. Chagin, A. Nowacki, and M. Kattan. Strategies for handling missing data in electronic health record derived data. EGEMS, 1(3):1035, 2013.
 [31] H. Xu, D. Luo, X. Chen, and L. Carin. Benefits from superposed Hawkes processes. In Proceedings of the 16th International Conference on Artificial Intelligence and Statistics (AISTATS), 2018.
 [32] H. Xu, D. Luo, and H. Zha. Learning Hawkes processes from short doublycensored event sequences. In ICML, 2017. arXiv:1702.07013.
 [33] Y. Yang, J. Etesami, N. He, and N. Kiyavash. Online learning for multivariate Hawkes processes. In NIPS, 2017.
 [34] K. Zhou, H. Zha, and L. Song. Learning social infectivity in sparse lowrank networks using multidimensional Hawkes processes. In Proceedings of the 16th International Conference on Artificial Intelligence and Statistics (AISTATS), 2013.
Supplementary Material
The supplementary material is organized as follows. In Section A we recall all the notation needed for the proofs. In Section B we restate the main statistical results from the paper and in Section C we provide proofs of these results. In Section D we state the optimization results and give proofs in Section E.
Appendix A Notation and Preliminary Remarks

negative log likelihood of complete data given th row .

: degree Taylor series approximation of .

negative log likelihood of missing data given th row . Loss function is unbiased in the sense that .

degree Taylor series approximation of



: fraction of data which is observed

:
Finally, we introduce additional notation which will be helpful in the proofs of Lemmas B.4 and B.5. First let denote the set of all monomials of degree . We represent an element as a list containing elements. An element in the list corresponds to the index of a term in the monomial (the list can potentially have repeated elements). For an example, the monomial can be represented as the list .
For a polynomial function we let denote the coefficient of the monomial in . Finally we define the order of a list to denote the number of unique elements in the list, so whereas .
Example
Consider the function . Using all the notation above, we can decompose as
where with corresponding coefficients , and .
Remark 1.
We next make several observations by applying the notation above to functions which appear in the likelihoods and . First, for a fixed we decompose the following function as a sum of monomials:
and note that
(A.1) 
We also have
and we similarly note that
(A.2) 
Next consider the function which appears in the missing data likelihood .
where . This observation will be important for our analysis because it allows us to leverage Equation A.1. Similarly we have
where allowing us to use Equation A.2.
Appendix B Statistical Results
We assume and . We take .
Theorem B.1 (Accuracy of ).
Suppose where . Then
for with probability at least .
The proof of Theorem B.1 relies on the following supplementary lemmas.
Lemma B.2.
Let . Then .
Lemma B.3 (Truncation error of ).
Suppose . Then
Lemma B.4 (Truncation Error of ).
Suppose . Then
Lemma B.5.
with probability at least .
Lemma B.6.
with probability at least .
Lemma B.7.
Let and . For any we have
and
with probability at least .
Appendix C Proofs of Statistical Results
c.1 Proof of Theorem b.1
Part 1: Controlling the Remainder
We set and let . Note that the loss functions are decomposable, i.e., . Since