1 Introduction
For more than forty years, social scientists have been developing datasets of events for the quantitative analysis of world politics. The last decade has witnessed a dramatic increase in activity in this area, much of it focused on automatic event detection for purposes of explaining and predicting political crises Beieler2016 . All of these efforts however, have used public information, such as newspaper or wire service reporting. Rather than directly measuring political activity, these systems can only count what reporters write about, which can vary over time and geography depending on many extraneous factors jenkins2016should . Together with the intrinsic challenges in automatic extraction, this has produced datasets that purport to track the same kind of events, such as political protests, but that are completely uncorrelated Hanna2010 . Moreover, some of the most important political activity is not immediately reported, and may not become publicly known until decades later, when formerly secret records are declassified. Even then, the sheer volume of these records can make it difficult even for the diligent researcher to identify individual events and assess their relative importance.
In this paper we study a new dataset of declassified documents and use statistical methods to identify and rank political events. Since 1973, the State Department has been using electronic records systems to preserve classified communications. The National Archives^{1}^{1}1Website: https://aad.archives.gov/aad/serieslist.jsp?cat=WR43 now makes available over 1.4 million declassified cables from 197377 and in addition, the metadata of more than 0.4 million other communications originally delivered by diplomatic pouch. They are all machinereadable and rich with metadata, creating many opportunities for statistical modeling.
(a) # of communications for TAGS CY  (b) # of communications for TAGS VS  
Count 


(c) # of communications for TAGS UNGA  (d) # of communications for TAGS FI  
Count 
Our goal is to explore modern statistical techniques that can automatically identify statistically interesting events in an important corpus of historical documents, which will continue to grow yearbyyear as millions more communications are declassified. For these communication streams, we are interested in studying and identifying “interesting” statistical patterns – we contend that these patterns correspond to heightened diplomatic activity; and validate our findings with standard reference works on the 1970s. A statistically interesting pattern can mean several things. Loosely speaking, this can correspond to sudden localized changes or abrupt “jumps” in communication traffic, regardless of the overall seriesspecific baseline activity (a communication stream may be very active or have very low traffic intensity overall). There can also be continuous periods in a communication stream, where the data lies consistently above a seriesspecific baseline that corresponds to a representative global activitylevel of that stream – these are “bursts” of activity in the temporal structure of the document streams that probably correspond with heightened diplomatic activity, such as the start or end of a war. An interesting event can also correspond to a heightened traffic intensity that plays out over longer periods, such as an increase over time, whether or not there are jumps.
When these communications were first entered in the State Department system, they were assigned one or more TAGS (Traffic Analysis by Geography and Subject), which indicate what countries or subjects each cable is related to. For example, “VS” signifies South Vietnam, and “SHUM” concerns human rights. By using these contentbased TAGS as the feature, we avoid the complication of language processing and focus on identifying statistically relevant activity patterns in the communication streams.
1.1 A brief exploratory description of the data
A glimpse of processed data in the form of communication streams is shown in Figure 1. The data shows that there is less traffic on weekends and holidays (including the end of the year). In addition, the number of communications sent in 1973 seems to be smaller compared to later years, due to fewer measurements. Since the overall (aggregated across all TAGS) number of communications sent across time did not have any systematic pattern indicative of events of historical importance, we decided to study the time series at a more granular level, by restricting to different types of TAGS. In Figure 1, panels (a)–(d) represent the communication streams, when restricted by TAGS type. Panels (a)–(c) show noticeable forms of increased activities in portions of the series – these are indicative of “interesting” historical events. For example, in panel (a) the increased activity in July 1974 corresponds to the Cyprus coup; in panel (b) the increase in number of diplomatic communications in April 1975 corresponds to the Fall of Saigon. Panel (c) shows multiple bursts in the number of cables, containing the particular TAGS UNGA (which stands for United Nations General Assembly) – interestingly, they all correspond to the annual United Nations General Assembly meetings. In addition to these visible bursts there seem to be some shorter periods of heightened activities, such as the smaller peaks for VS (South Vietnam) a year after the fall of Saigon corresponding to the ensuing refugee crisis.
In contrast to panels (a)–(c) in Figure 1, panel (d), for cables related to Finland (FI), does not seem to show any period of heightened activity during the time period under consideration. These prototypes are representative of the different TAGSspecific series: Exploratory analyses of the database of TAGS specific communication streams suggest that there are several series with some “interesting event” (as in Panels (a)–(c)), while others seem to be less active (as in panel (d)). A first goal of our work is to quantitatively define traits that separate communication streams like the figure in panel (d) from those in panels (a)–(c). We develop statistical methods that can mine these (TAGS specific) timeseries and identify communication streams that exhibit statistically interesting activities in them. Once we identify these interesting communication streams, we develop algorithms that perform a deeper investigation of each series and identify time intervals where the signal undergoes abrupt localized changes in communication traffic. The informal ideas described above are made more precise in Section 2 of this paper.
Exploratory analysis suggests that changes in the proportion of a particular TAGS appearing in a communication stream are better representatives of identifying whether a period is active or not, as compared to tracking the corresponding counts. Due to the noticeable difference in the number of cables that were communicated over the weekdays and lowtraffic days, the communication streams shown in Figure 1 seem to be a superposition of a high traffic series and a lowtraffic (weekend and holiday) series. As a preprocessing step, we filtered out the days where the total number of cables being communicated were very small – they lead to more reliable estimates of proportions.
(a)TAGS FI, pvalue = 0.181  (b)TAGS OSCI, pvalue = 0.005  
Proportion 


(c) TAGS PFOR, pvalue = 0  (d) TAGS PINT, pvalue = 0  
Proportion 
2 Statistical Methodology
For the reader’s convenience, we first present an outline of the main statistical approaches that we discuss in this paper. Section 2.1 addresses our first question: how do we determine whether a communication stream, among several hundreds, is interesting or not? We address this as a statistical hypothesis testing problem, where we consider the problem of whether a communication stream is generated from a homogeneous binomial process – in other words: is the proportion of documents containing the particular TAGS uniform across time? We select a collection of TAGS for which there seems to be some form of a statistically interesting event, and explore their geometry in further detail. If denotes the proportion of cables containing a TAGS at time , we estimate the function with a piecewise constant signal – this is performed by using a regularized negative loglikelihood criterion based on the fused Lasso penalty TSRZ2005 or its counterpart (boysen2009consistencies, ; Johnson2013fldp, )
. For efficient computation, we develop new fast algorithms for these penalties by adapting existing work for the least squares loss function to the case of the generalized linear model likelihood, studied herein – see Section
2.2. The piecewise constant segments lead to localized changes in the underlying signal – we use hypothesis testing ideas based on sample splitting wasserman2009high to rankorder the strengths of the jumps based on the associated pvalues – see Section 2.3. The locations of discontinuities or jumps of the signal are aggregated to obtain locally contiguous subintervals of heightened diplomatic communications, which we call bursts, adopting the terminology coined by kleinberg2003bursty in a different context (Section 2.4). Finally in Section 3, we discuss preliminary results on estimating the underlying signal with more flexible local models, beyond piecewise constant segments.2.1 Identifying Interesting Communication Streams
We develop a framework to identify communication streams that exhibit some form of a heightened or in other words, statistically interesting activity. Consider a TAGSspecific series , where, denotes the number of documents containing the specific TAGS among cables, with proportion . We will assume that the conditional distributions of ’s are all independent across . We ask:
Is there any evidence of (localized) heightened intensity of the proportions, compared to a baseline model, where all proportions are the same?
To measure a localized change (increase) in intensity, we fix a window of size and consider all the points in the neighborhood of a time point , given by: (we pick days in our experiment). The average proportion in this neighborhood: , is a measure of communication traffic around the reference time point . We say that a large value of compared to a baseline value (for example, the global proportion), indicates the presence of an intense activity of some form; and declare such an occurrence to be statistically interesting. We hypothesize such an event to be associated with an event of historical interest; and subsequently validate this belief by factoring in the insights of a historian or social scientist.
Formally, we consider a simple testing problem with : versus there exists an such that is larger than the (global) average proportion. Inspired by popularly used scan statistics glaz2001scan
, we propose to use the following test statistic:
(1) 
where,
is the (estimated) global proportion of the signal estimated under the null hypothesis;
is an estimate of ; andis the standard deviation of
evaluated under the null (). measures the strength of a locally contiguous period of heightened activity–the larger this value, the more pronounced is the localized traffic. We use a permutation based approach to compute the distribution of under the null. Figure 2 shows different communication streams with their associated pvalues computed using the framework described above; and varying levels of activity in the different representative communication streams – a large pvalue for panel (a) (corresponding to TAGS FI) signifies a lack of interesting activity in this series – this aligns with an expert’s understanding of historical events. A table providing a summary of how many cables survive different pvalue thresholds are provided in Table 1 (Appenfix).We also used a teststatistic depending upon the likelihood ratio test: we replaced (in display (1
)) by the local likelihood ratio derived from the binomial distribution; assuming that
are independently distributed for . The results obtained via this likelihood based approach was quite similar to that obtained from the modelfree testing procedure described above; and are hence omitted for the sake of brevity.2.2 Identifying Jumps in Communication Streams
The method in Section 2.1 simply identifies whether something statistically interesting, in the form of a heightened localized activity (for example), is occurring somewhere in a communication stream. It does not inform us about the structure or precise location of such an activity. We proceed to explore some finer structural properties of the series. Inspired by popularly used signal segmentation/estimation methods TSRZ2005 ; Johnson2013fldp ; mammen1997locally in statistical signal processing, we seek to identify breaks or jumps in a piecewise constant approximation of the signal – these changes are more localized than the (global) deviations studied in Section 2.1.
We use the notation of Section 2.1. We assume herein that , where, denotes the probability of success. Our model assumes that the parameters ’s are piecewise constant. Locations where the underlying signal exhibits a discontinuity will be called a “jump” in the communication stream. Instead of working directly with the proportion ’s we will use the logisticlink function: and will treat as a natural parameter. This leads to the following regularized estimation problem:
(2) 
where, , the negative logarithm of the likelihood is the datafidelity term and is the regularizer. encourages the estimated ’s (and hence the proportion ’s) to be piecewise constant and the regularization parameter controls the amount of shrinkage. Two common examples of that we use herein are TSRZ2005 ; mammen1997locally ; boysen2009consistencies ; Johnson2013fldp :

segmentation (Fused Lasso): , which penalizes the total variation of a signal, which may also be thought as a softversion of the number of the number of jumps in .

segmentation: Here we take , which penalizes the number of jumps in the signal .
We assume above and in the discussion below, that the time points are equally spaced. If they are not equispaced, the penalty function needs to be adjusted in a straightforward fashion as discussed in Section A.1.
For the penalty function , Problem (2) is a convex optimization problem. The penalty on the successive differences in , is commonly referred to as the fused lasso or totalvariation penalty TSRZ2005 ; mammen1997locally . This seminorm induces shrinkage along with sparsity on the coefficient differences . However, the based penalty often overestimates the number of jumps when the tuning parameter is chosen so as to obtain a model with good datafidelity – this is especially true if the tuning parameter is chosen based on a held out test set to minimize test error. This is due to the shrinkage effect of the penalty, which severely penalizes large values of the jumps . To obtain a model with fewer jumps, the regularization parameter needs to be made larger – this however, may lead to a model where some of the important jumps are missed. Many of these limitations can be overcome by using a based penalty boysen2009consistencies which directly penalizes the number of jumps and is less agnostic to the precise value of the jump. The caveat however, is that the resulting optimization problem becomes nonconvex and discrete optimization methods are required to obtain the global minimum of such problems boysen2009consistencies ; Johnson2013fldp . Developing global optimization algorithms for Problem (2) (for the penalty ) along the lines of mazumder2017discrete ; bertsimas2016best is beyond the scope of the current paper. Herein, we describe fast and robust algorithms to obtain good quality solutions to Problem (2). Our proposal is motivated by firstorder optimization based algorithms nesterov2004introductorynew pioneered in the convex optimization community. Loosely speaking, these are iterative methods that can be used to obtain high quality approximate solutions for convex optimization tasks, compared to offtheshelf interior point methods that are difficult to scale to large problems. These methods apply to both and – however, there are certain subtleties as we describe below.
2.2.1 Model Fitting: Optimization Algorithm
Problem (2) is of the composite form, i.e., the objective function can be written as the sum of a smooth convex function and a nonsmooth penalty function . We will apply proximal gradient descent methods fista09 for this problem. The negative loglikelihood function is continuously differentiable and satisfies:
(3) 
with . This follows by noting that the th coordinate of is: and is a diagonal matrix with the th diagonal entry being:
(4) 
Hence, the largest eigenvalue of
, i.e., , which justifies the choice of , as above. For a fixed , our algorithm performs the following updates (for all ):(5) 
This leads to a decreasing sequence of objective values for . We now study the fate of the sequence , depending upon the choice of .
The fused lasso penalty ()
We first consider the case where the regularizer and hence the function is convex in . The sequence converges to a minimum of Problem (2) with the penalty function . More precisely, it can be shown that fista09
(6) 
where, is an optimal solution to Problem (2). Display (6) states that the sequence converges to a minimizer of Problem (2) with a rate of . In fact, this rate can be improved further under a minor additional assumption. Let , i.e., the smallest diagonal entry of the Hessian of . As long as ’s are uniformly bounded and then . The rate of convergence is linear; and is given by:
(7) 
where,
(8) 
Note that subproblem (5) requires solving a problem of the form:
(9) 
which can be solved very efficiently via Dynamic Programming Johnson2013fldp with cost for – one can solve instances with a few thousands in a few milliseconds.
The segmentation penalty ()
The algorithm above, also applies to the regularizer . In update (5) we set to . This requires solving
(10) 
which can also be computed efficiently using the dynamic programming algorithm of Johnson2013fldp . Describing the properties of this sequence is subtle since the associated optimization problem (2) is nonconvex. Following bertsimas2016best it can be shown that the sequence is decreasing, bounded below^{2}^{2}2We will assume that the function is bounded below and the minimum exists. and it converges to (say), which may not be a global minimum. We say that is a firstorder stationary point for Problem (2) if it satisfies the following fixed point equation:
is said to be an accurate first order stationary point, if . Following the convergence analysis in bertsimas2016best [Theorem 3.1], we obtain the following finitetime convergence rate of to a first order stationary point:
(11) 
Note that this algorithm may not reach the global minimum of the version of Problem (2). However, in practice, it reaches a high quality solution quite fast.
2.2.2 Estimated Signal
To gather some intuition about the behavior of the segmentation methods described above, we consider a synthetic example in Figure 3 – here the underlying (true) signal is piecewise constant with three levels up to time point , there is a right discontinuity at after which it becomes linear ^{3}^{3}3More specifically, data is generated by , where for ; for ; for ; for .. The figure presents the signal estimates (for both the and penalties) at the crossvalidated choices of the tuning parameter – we use fold cross validation FHT09new which is also used in the R package genlasso (Since we want to ensure each fold is representative of the time series, instead of randomly assigning points to a fold, we systematically assign points by placing every th point into the same fold). For both segmentation schemes, the estimated signals serve as good (overall) approximations of – however, there are some subtle differences. First of all, the segmentation scheme leads to biased estimates and the bias becomes quite prominent in estimating the jump at the centre of the signal. This behavior is not present for the scheme. In addition, the estimates for the linear component (at the right) also differ across the and schemes. The regularizer leads to a fewer number of segments (here three) compared to the penalty which has several smaller jumps.
regularization  
Proportion 


regularization  
Proportion 
Figures 4 and 5 show the estimated signal proportions , where the ’s are obtained from Problem (2). Both the penalty functions do a good job in estimating a piecewise constant version of the underlying signal – the scheme leads to fewer jumps than its counterpart, for a comparable datafidelity. The figures also show fitted signals for a few other values of around the crossvalidated choice at the center (
increases as one moves down the rows): we include the tuning parameter selected by the onestandard error rule
FHT09new (see also the R package genlasso). We can see that as decreases, the algorithm captures a more granular structure of the data and estimates more jumps.segmentation  segmentation  

Proportion 

Proportion 

Proportion 
segmentation  segmentation  

Proportion 

Proportion 

Proportion 
2.3 A deeper investigation of Jumps
2.3.1 How intense is a Jump?
The procedure in Section 2.2 provides an estimated piecewise constant signal . In particular, this can be used to identify locations where the signal changes: . Of course, this jump depends upon the choice of the tuning parameter, the penalty used and also the underlying signal. A jump estimated by the or segmentation procedure may reflect (a) a discontinuity in the signal – in this case, the signal is well approximated by locally constant segments with pieces adapting to the data (b) a localized trend in the signal, as we saw in Figure 3 – a jump here is a consequence of the slope of and not a discontinuity in the signal .
Given an estimate of a scholar accustomed to analyzing events through a close reading of historical documents may ask:

Which of these jumps might be important or are indicative of a historical event of interest?

Can one obtain a rank ordering of the jumps based on their intensities?
We formalize this question as follows: given an estimate of and a set of candidate jumps, can we identify jumps that are strong enough? Ideally, we would like a simple measure that associates a score to the strength and size of a jump selected by the estimation procedure – this would help us select a smaller set of jumps that merit closer scrutiny. Towards this end, we use a sample splitting^{4}^{4}4Due to the large number of samples, sample splitting does not significantly reduce the size of the training dataset. procedure wasserman2009high : a subsample of size 50% of the data is used for estimating the location of the jumps and the remaining held out part of the data is used to associate a pvalue score (the method is described below) to each jump identified in the first stage. This method is simple, intuitive and provides a natural method to rank the jumps in a communication stream. Using this scheme, one can potentially reduce the number of jumps selected by the fitting procedure by screening out jumps with pvalues larger than a userdefined prespecified threshold.
We describe our approach with reference to the penalization procedure, though the idea will also apply for the penalized estimate. Suppose a candidate location for the change point is estimated based on the first part of the sample (used for estimating the signal). The test statistic is evaluated on the held out sample. We take a neighborhood of size centered at , and denote the time points on the left of as and those on the right of as . Let us assume that for are all equal to ; and for are all equal to . We then test the null hypothesis () that the proportions on the left and right parts of are equal: ; versus the alternative () that . We use the likelihood ratio test statistic for this purpose. To compute the null distribution, we used a two step procedure. We first identified segments of the time series which did not overlap with any candidate change point location (i.e., parts of the series where the estimated signal was constant for stretches of size at least ) of the time series. Based on the regions thus selected (i.e., the locally constant stretches of the signal), we simulated the null distribution of the test statistic by using a permutation test. These 2sided pvalues were consequently used as a measure of intensity of every jump.
Note that a candidate jump estimated by the signal estimation procedure at the crossvalidated choice of the tuning parameter need not necessarily correspond to a jump with a low pvalue. A small pvalue indicates that the intervals to the left and right of have different proportions^{5}^{5}5A jump obtained from the estimated may be due to a linear rise in the signal which need not correspond to a significant change in local proportions. Our experiments indicate that jumps in that correspond to gradual linear rises in the signal, have higher pvalues associated with them when compared to sudden or abrupt changes in . Thus the pvalues can be used to (a) devise a scoring mechanism to rank order multiple jumps observed in a series and/or (b) prune out redundant jumps and identify ones that exhibit a strong difference in proportions between the left and right intervals. Figure 6 shows the communication stream for TAGS CVIS and the estimated signal obtained via segmentation. We also computed the pvalues for each potential jump location as suggested via the segmentation fit. We pruned them and refitted the model based on their pvalue scores. It helps to interpret these patterns alongside Figure 9 (TAGS CVIS) which presents more flexible fits of the underlying signals for this communication stream. Figure 7 shows additional examples interpreting the pvalues associated with the different jumps. Figure 7 suggests that the pvalues are indicative of whether a jump is due to a shift in the piecewise constant level or a linear trend – the pvalues are larger when there is a linear trend rather than a sharp jump (as in a piecewise constant signal). Figure 9 shows a more flexible approximation of this communication stream which provides further insights into the patterns of TAGS US series. To validate the intuition gathered above, we consider a synthetic example in Figure 8 (with the same data as in Figure 3) – here we observe that the pvalues tend to be larger for jumps in the right part of the signal – these jumps in the piecewise constant segments result from estimating a linear trend (that appears at the right of the series) with piecewise constant segments. Note that the pvalues associated with the first three jumps (at the left of the signal) are quite small – they correspond to jumps in the underlying piecewise constant signal.
In passing we note that it is also possible to perform multiple testing procedures lehmann2006testing to attach error rates to a family of jumps. In particular, we can use Family Wise Error Rates or False Discovery Rate (FDR) control procedures to return a collection of candidate jumps with a certain prescribed control on the error rate lehmann2006testing . However, our goal here is to use the sample splitting and pvalue framework to perform an exploratory analysis of the strengths of the different jumps – we do not pursue an indepth study of multiple testing in this work.
No pruning  Pruning with  
Proportion 


Pruning with  Pruning with  
Proportion 
TAGS ENRG  TAGS US  

Proportion 
Proportion 
2.4 From Jumps to Bursts
Section 2.2 presents methods to detect several jumps in a signal, and also presents an indepth investigation of how to associate scores (pvalues) to each of the detected jumps using a sample splitting scheme. We now explore if it is possible to summarize a single communication stream (corresponding to a specific TAGS) with a score – a measure that should ideally contain in it information pertaining to the strengths and stretches of the different jumps. Towards this end, following the terminology introduced by Kleinberg kleinberg2003bursty we formalize the notion of a “burst”. The general approach pursued in this paper and the models used differ from kleinberg2003bursty . Informally speaking, a burst corresponds to a stretch of time where a communication stream depicts traffic which is larger than a baseline value. We present a computational scheme to estimate such bursts for a TAGS specific communication stream.
2.4.1 Computation of the strength of a Burst
As a starting point, we consider an estimate of a baseline proportion (we discuss how to compute this below) that is specific to a communication stream. A “burstiness period” or simply burst corresponds to a time interval where the estimated signal lies above the baseline value and is given by , where . Following kleinberg2003bursty , we define the strength of the burst as the logarithm of the likelihood ratio (here, the numerator is the likelihood of the signal and the denominator is that evaluated at the baseline) given by: where denotes the likelihood at time . Note that the baseline is specific to a communication stream and the score represents a deviation from this global baseline. is different than the magnitude of a jump given by – it takes into account the deviation of from the baseline as well as the duration of the burst given by the length of . A large value of means that a large part of the likelihood is explained by deviations from the baseline, and therefore, corresponds to a strong burst. Note that each TAGS can have multiple bursts leading to multiple intervals – each with an assigned strength .
Choice of baseline: The baseline value should be representative of the behavior of the TAGspecific communication stream. The global proportion of a communication stream is a reasonable choice. We set to be one standard deviation larger than the global proportion
A robust estimate like the median can also be used instead of the average. In our experiments we found that the topranked slots were relatively agnostic to the choice of the baseline .
2.4.2 Interpretation of Bursts
Table 2 presents the top thirty bursts, with the start and end dates, and the date with the highest value. A close study of the content of the cables shows that not all of these bursts correspond with what scholars would recognize as an even of historical importance. After all, the cable TAGS that diplomats used do not necessarily correspond with diplomatic activity. For instance, the second biggest burst is made up of cables related to transportation (ETRN) a TAGS that was commonly used, and overused, from when we begin to have records continuing until 1974, when diplomats’ use of this TAGS was largely discontinued. The biggest burst, for CVIS (visas), has a similar pattern (as shown in Figure 6). But in this case, it appears to reflect a decision by archivists to stop preserving records related to visas langbart2007 . To the model, both of these look like bursts, but they simply reflect administrative procedures rather than historical events.
The bursts that follow, on the other hand, appear to correspond well to historical events. The next ten include the Carter administration’s prioritization of human rights (SHUM), Anwar Sadat’s surprise visit to Israel (PGOV), the Southeast Asian “Boat People” crisis (SREF), the U.S. withdrawal from the International Labor Organization (PORG), the conclusion of the Panama Canal Treaty (PDIP), the 1973 Yom Kippur War (XF, for Middle East), Portugal’s withdrawal from Angola (AO), and the 1974 crisis over Cyprus (CY). All are included in each of the four standard reference works we consulted brune2003chronological ; flanders1993dictionary ; jentleson1997encyclopedia ; deconde2002encyclopedia .
A systematic evaluation of hundreds of bursts for historical significance lies outside the scope of this paper. But the relative proportion of recognized historical events appears to diminish as one examines smaller bursts, like the ones ranked in the range 1322. They include the denouement of the Vietnamese War (VM and VS), the OPEC oil embargo (ENRG), the Vladivostok summit (OVIP), and negotiations to end white rule in Rhodesia (RH). But there are also largely unrecognized events, like a 1975 UN General Assembly debate over the command of foreign military forces in South Korea, that would appear to merit closer scrutiny. The identification of such unstudied episodes, no less than rankordering wellknown events, is valuable for historical scholarship.
3 Beyond Piecewise Constant Segments
A major focus of Section 2 was on approximating a communication stream with a piecewise constant signal. This framework does help us answer some of the major datadriven questions of interest to a political scientist, based on a first order approximation of the communication streams. We now investigate more flexible signal approximations that provide us insights into the finer behavior of the signals. A natural extension of a piecewise constant estimate is a piecewise linear estimate. However, there are subtleties in incorporating this structure into our framework, as we discuss below.
To settle ideas, let us consider the usual signal denoising problem with data: , for where, . We seek to estimate such that it is piecewise linear. In this vein, it is common to use the trendfiltering approach kim2009ell_1 ; tibshirani2014adaptive with regularizer to obtain a signal with piecewise linear segments:
The penalty function encodes the norm on the discrete second order derivative of the signal assuming that the time points are all equally spaced. can be interpreted as a convexification of its version: that counts the number of different piecewise linear segments.
Our situation is different from the denoising example outlined above. Since we are working under the modeling assumption: with , imposing a trend filtering penalty on so as to maintain piecewise linearity will lead to a nonconvex optimization problem due to the nonlinear dependence of on . Thus, instead of enforcing the sequence to be piecewise linear, we allow the latent parameters to be piecewise linear – this leads to a computationally tractable estimation framework. Encouraging to be piecewise linear enables to be more flexible than a piecewise constant signal. Towards this end, we propose a simple adaption of the estimation criterion in Problem (2) by setting the regularizer . Figure 9 shows the results of estimates obtained from some communication streams using the trend filtering penalty. If the time points are not equally spaced, then this penalty can be modified appropriately – see for example kim2009ell_1 and also Section A.1.
Computation
The proximal gradientstylized algorithm update (5) can be adapted to the setting described above with . We use the specialized interior point solver^{6}^{6}6We use the Rpackage wrapper available from https://github.com/hadley/l1tf of kim2009ell_1 – this works quite nicely for the problem sizes encountered in this paper. The resulting Problem (2) is convex and the sequence (5) leads to the optimum of the problem. The convergence rates outlined in (6) and (7) will also apply to this problem. If we set , the resulting Problem (2) becomes a challenging nonconvex optimization problem – in this case, there is no analogue of the highly efficient dynamic programming implementation that was available for . Thus, for the case where we desire to be piecewise linear, we confine our study to the choice of the convex trend filtering regularizer.
Acknowledgements
R. Mazumder was supported by ONR (grant # N000141512342) and NSFIIS (grant # 1718258) for support. The authors will like to thank: David Blei, David Madigan, Shawn Simpson for helpful comments and encouragement. The authors will also like to thank the workshop participants of “Famine and Feast–International Historical Research in the Digital Age” (London, UK; 2015) for comments on the work. Preliminary results from this article appeared in a media article on buzzfeed.com ^{7}^{7}7Article www.buzzfeed.com/josephbernstein/canacomputeralgorithmdothejobofahistorian?.
CY  ENRG  
Proportion 


VS  SHUM  
Proportion 
Appendix A Appendix
Significance level  number of TAGS for test (1) 

0.1  914 
0.01  768 
0.001  622 
0.0001  509 
0.00001  391 
a.1 Problem (5) with irregularly spaced time points
If the time points are irregularly spaced then the penalty functions need to be adjusted accordingly. The fused lasso penalty function becomes: where, denotes the time difference between time point and the next time point, indexed by . For this choice, the associated proximal map i.e., Problem (5) needs to be modified to:
(12) 
where, . More generally, if we consider the trend filtering example with varying time intervals, we get an instance of Problem (12) with
We use the Alternating Direction Method of Multipliers (ADMM) procedure boydadmm1 which performs the following decomposition: and obtains the Augmented Lagrangian:
for some choice of . The usual ADMM approach performs the following sequence of updates:
(13)  
where, in the update wrt other variables remain fixed, and the same applies for the update wrt . We refer the reader to boydadmm1 for details pertaining to the convergence of this algorithm and choices of . We note that the update wrt in display (13) can be solved quite easily via solving a system of linear equations:
Note that is a bidiagonal matrix when corresponds to the weighted fused lasso penalty and a tridiagonal matrix when it corresponds to the weighted trend filtering penalty. The inverses in each of these cases can be computed with cost and (respectively) BV2004 ; kim2009ell_1 – furthermore the inverse can be computed once (at the onset) as the matrix does not change across iterations. The update wrt in display (13) requires a solving the following problem:
where, . A solution to the above problem is given by the familiar softthresholding FHT09new operation where, . The sequence of updates in (13) are performed till some form of convergence criterion is met boydadmm1 .
TAGS  meaning  start  end  peak  Burst Strength  

1  ETRN  Economic AffairsTransportation  19730702  19740809  19730928  5146.05 
2  CVIS  Consular AffairsVisas  19730702  19750102  19740628  4839.35 
3  SHUM  Social AffairsHuman Rights  19770119  19771230  19771118  2872.02 
4  US  United States  19760128  19770916  19760415  2516.03 
5  PGOV  Political AffairsGovernment  19770603  19771230  19771118  2484.57 
6  SREF  Social AffairsRefugees  19750422  19760720  19760602  1662.58 
7  SOPN  Social AffairsPublic Opinion and Information  19761126  19771230  19770826  1597.14 
8  PORG  Political AffairsPolicy Relations With International Organizations  19770615  19771230  19771111  1547.35 
9  PDIP  Political AffairsDiplomatic and Consular Representation  19770524  19771230  19770902  1462.93 
10  XF  Middle East  19731009  19731219  19731016  1453.76 
11  AO  Angola  19751108  19760223  19751110  1439.58 
12  CY  Cyprus  19740715  19740729  19740720  1378.79 
13  VM  Vietnam  19771011  19771230  19771012  1365.45 
14  PDEV  Political AffairsNational Development  19770613  19771230  19770831  1344.70 
15  VS  Vietnam (South)  19730702  19750606  19750425  1150.46 
16  UNGA  UN General Assembly  19750819  19751213  19751107  1044.29 
17  CARR  Consular AffairsAmericans Arrested Abroad  19770601  19771230  19770628  951.98 
18  MCAP  Political AffairsMilitary Capabilities  19730702  19740815  19740703  903.26 
19  ENRG  Economic AffairsEnergy  19731108  19740221  19740125  760.64 
20  PBOR  Political AffairsBoundary and Sovereignity Claims  19770701  19771230  19771109  685.75 
21  OVIP  OperationsVIP Travel Arrangements  19741009  19741109  19741031  607.17 
22  RH  Rhodesia  19760901  19771230  19770831  569.89 
23  AEMR  AdministrationEmergency and Evacuation  19750328  19750512  19750428  524.59 
24  MPLA  Popular Movement for the Liberation of Angola  19751107  19760224  19760218  507.98 
25  MSG  Marine Security Guards  19760902  19771230  19771128  507.27 
26  OREP  OperationsCongressional Travel  19761027  19761118  19761102  481.08 
27  PRG  Provisional Revolutionary Government of South Vietnam  19750116  19750206  19750203  470.51 
28  MNUC  Military and Defense AffairsMilitary Nuclear Applications  19770311  19771230  19770822  421.53 
29  UNGA  UN General Assembly  19740905  19741205  19741010  417.20 
30  CB  Cambodia (Khmer Republic)  19730702  19750521  19750416  370.14 
SA  TQ  FI  

Proportion 


SW  AR  OSCI  

Proportion 

PFOR  SREF  OEXC  
; rapid change 
Proportion 

PINT  OCON  SF  
; gradual change 
Proportion 
References
 [1] Amir Beck and Marc Teboulle. A fast iterative shrinkagethresholding algorithm for linear inverse problems. SIAM J. Imaging Sciences, 2(1):183–202, 2009.
 [2] John Beieler, Patrick T. Brandt, Andrew Halterman, Philip A. Schrodt, and Erin M. Simpson. Generating political event data in near real time: Opportunities and challenges. In R. Michael Alvarez, editor, Computational Social Science, pages 98–120. Cambridge University Press, 2016.
 [3] Dimitris Bertsimas, Angela King, and Rahul Mazumder. Best subset selection via a modern optimization lens. Annals of Statistics, 44 (2):23–42, 2016.
 [4] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. Number 3(1). Now Publishers, 2011.
 [5] Stephen Boyd and Lieven Vandenberghe. Convex Optimization. Cambridge University Press, Cambridge, 2004.
 [6] Leif Boysen, Angela Kempe, Volkmar Liebscher, Axel Munk, and Olaf Wittich. Consistencies and rates of convergence of jumppenalized least squares estimators. The Annals of Statistics, pages 157–183, 2009.
 [7] Lester H Brune and Richard D Burns. Chronological History of U.S. Foreign Relations: 19331988. Routledge, 2003.
 [8] Alexander De Conde, Richard D Burns, Fredrik Logevall, and Louise B. Ketz. Encyclopedia of American foreign policy. Scriber, 2002.
 [9] Stephen A Flanders and Carl N Flanders. Dictionary of American foreign affairs. Macmillan Library Reference, 1993.
 [10] Joseph Glaz, Joseph I Naus, Sylvan Wallenstein, Sylvan Wallenstein, and Joseph I Naus. Scan statistics. Springer, 2001.
 [11] Alex Hanna. Assessing gdelt with handcoded protest data. www.badhessian.org, 2014. (accessed: July 29, 2016).
 [12] Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The Elements of Statistical Learning, Second Edition: Data Mining, Inference, and Prediction (Springer Series in Statistics). Springer New York, 2 edition, 2009.
 [13] J Craig Jenkins and Thomas V Maher. What should we do about source selection in event data? challenges, progress, and possible solutions. International Journal of Sociology, 46(1):42–57, 2016.
 [14] Bruce W Jentleson, Thomas G Paterson, and Nicholas X Rizopoulos. Encyclopedia of US foreign relations, volume 2. Oxford University Press, USA, 1997.
 [15] Nicholas A. Johnson. A dynamic programming algorithm for the fused lasso and l0segmentation. Journal of Computational and Graphical Statistics, 22(2):246–260, 2013.
 [16] SeungJean Kim, Kwangmoo Koh, Stephen Boyd, and Dimitry Gorinevsky. ell_1 trend filtering. SIAM review, 51(2):339–360, 2009.
 [17] Jon Kleinberg. Bursty and hierarchical structure in streams. Data Mining and Knowledge Discovery, 7(4):373–397, 2003.
 [18] David Langbart, William Fischer, and Lisa Roberson. Appraisal of records covered by N159073P. National Archives, 2007.
 [19] Erich L Lehmann and Joseph P Romano. Testing statistical hypotheses. Springer Science & Business Media, 2006.
 [20] Enno Mammen, Sara van de Geer, et al. Locally adaptive regression splines. The Annals of Statistics, 25(1):387–413, 1997.
 [21] Rahul Mazumder and Peter Radchenko. The discrete dantzig selector: Estimating sparse linear models via mixed integer linear optimization. IEEE Transactions on Information Theory, 2017.
 [22] Yurii Nesterov. Introductory Lectures on Convex Optimization: A Basic Course. Kluwer, Norwell, 2004.
 [23] R. Tibshirani, M. Saunders, S. Rosset, J. Zhu, and K. Knight. Sparsity and smoothness via the fused lasso. Journal of the Royal Statistical Society, Series B, 67:91–108, 2005.
 [24] Ryan J Tibshirani. Adaptive piecewise polynomial estimation via trend filtering. The Annals of Statistics, 42(1):285–323, 2014.
 [25] Larry Wasserman and Kathryn Roeder. High dimensional variable selection. Annals of statistics, 37(5A):2178, 2009.
Comments
There are no comments yet.