text:
A Grammar-based Compressed Representation of 3D Trajectories

Much research has been published about trajectory management on the ground or at the sea, but compression or indexing of flight trajectories have usually been less explored. However, air traffic management is a challenge because airspace is becoming more and more congested, and large flight data collections must be preserved and exploited for varied purposes. This paper proposes 3DGraCT, a new method for representing these flight trajectories. It extends the GraCT compact data structure to cope with a third dimension (altitude), while retaining its space/time complexities. 3DGraCT improves space requirements of traditional spatio-temporal data structures by two orders of magnitude, being competitive for the considered types of queries, even leading the comparison for a particular one.

Introduction
Geopositioned data is ubiquitously and continuously generated to describe different types of trajectories; e.g. routes of professional transportation vehicles or our daily running paths. Obviously, large and varied trajectory datasets are being consolidated, and they are exploited for different and innovative purposes. Disregarding their final application, managing trajectory datasets poses many challenges that have attracted much research efforts.
A prominent domain that demands efficient trajectory management is Air Traffic Management (ATM). ATM systems analyze very large flight-related datasets to make decisions to improve air traffic performance, reducing costs, or making safer and environmentally friendly airspaces. Currently, ATM services are evolving to support and leverage "next generation" technologies like Automatic Dependent Surveillance-Broadcast (ADS-B). ADS-B is a surveillance technology in which aircrafts determine flight parameters (latitude, longitude, altitude, etc) via navigation systems, and broadcast them to ground stations, that then deliver this data to ADS-B providers; e.g. the OpenSky Network [16], that is the provider of the ADS-B datasets used in our experiments.
ADS-B has been progressively adopted by many aircraft manufacturers, and more ground stations have been deployed around the world. It has increased ADS-B coverage, and also the size of ADS-B datasets, whose storage and querying has become more difficult. Storage issues were first addressed using columnar compression [20,22]. Although their numbers are moderately successful, the resulting representations can not be efficiently queried. More recently, a compressed index for ADS-B (called ADS-BI) has been proposed [21]. It performs block partitioning and stores descriptive metadata about the block to enable some types of queries. Block contents are then encoded by columns using universal compression (e.g. gzip or p7zip), reporting competitive numbers. Although ADS-BI resolves some type of queries by time or 2D-position, it does not support altitude-based searches, which is highly desirable for ATM systems; for instance, when a controller looks for aircrafts flying at certain flight level in a given region.
Therefore, our main objective is to propose a data structure that allows 3D trajectories to be effectively compressed, and searches to be performed by time and/or any of the three positional dimensions. It is not a new problem [7], and some researches have been previously published about 2D (latitude, longitude), and 3D (including altitude) trajectory management. Data structures like 3DRtree [19], HR-tree [13], the MVR-tree [17], or PIST [2] have been successfully used for many years, but currently show scalability issues when they are used to manage larger trajectory datasets. The Douglas-Peucker algorithm [8] has been used to make trajectories more compact; other examples are dead reckoning [18], TrajStore [6] and Trajic [15].
Our approach, called 3DGraCT, proposes a new compact data structure that stores and indexes 3D trajectories in compressed space. 3DGraCT enhances GraCT [14] to manage altitude information, and also to enable query resolution by this dimension. Our experiments, using different-size ADS-B datasets, show that 3DGraCT improves space requirements of traditional spatiotemporal data structures by two orders of magnitude, and competes with them in query performance, leading the comparison for queries asking for large time intervals.
2 Background k 2 -tree. The k 2 -tree [5] is conceptually an unbalanced k 2 -ary tree constructed from a binary matrix by recursively subdividing the matrix into k 2 submatrices of the same size, if k = 2, it is a space/time efficient version of a region quadtree. First, the original matrix is divided into k 2 submatrices of size n 2 /k 2 , being n×n the size of the matrix. Each of these submatrices generates a child of the root node whose value is 1, if there is at least one 1 in the cells of that submatrix, and 0 otherwise. The subdivision continues recursively for each child with value 1 until a submatrix full of 0s is found or the cells of the original matrix (i.e., submatrices of size 1 × 1) ar reached. Figure 1 shows an example of this subdivision (left) and the resulting conceptual k 2 -ary tree (right up) for k = 2.
The k 2 -tree is stored using two bitmaps T and L (see Figure 1). T stores all the bits of the k 2 -tree, except those in the last level, following a level-wise  traversal: first the k 2 binary values of the children of the root node, then the values of the second level, and so on. L stores the last level of the tree. k 3 -tree. The k 2 -tree can be generalized to deal with a three-dimensional binary cube, instead of a two-dimensional binary matrix. It can be trivially done by extending the space partitioning, while maintaining the representation techniques used for k 2 -trees. Thus, each 1 in the binary cube of the k 3 -tree [1] represents a tuple x, y, z , where (x, y) are the coordinates in the 2D space, and z is the altitude. It is possible to obtain efficiently the value of a cell, a cube, or slices of the cube, by just performing rank and select operations [10] over T and L.
Re-Pair. Re-Pair [12] compresses a sequence by recursively substituting pairs of symbols by a new one. Given a sequence of integers I (called terminals) the compression process is as follows: (1) it obtains the most frequent pair of integers ab in I; (2) it adds rule W → ab to dictionary R, where W is a new symbol not present in I (called a non-terminal); (3) every occurrence of ab in I is replaced by W , and (4) it repeats steps 1-3 until all pairs in I appear only once (see Figure 2). The resulting sequence after compressing I is called C.
GraCT. GraCT [4] is a compact data structure to represent and query trajectories of moving objects in a free space of two dimensions. It requires that all objects declare their positions at regular time instants (e.g. each minute), but interpolation is used when an object does not inform its position in a given instant. GraCT uses a raster model to represent the space; i.e. it is divided into cells (squares) of a fixed size, and it is assumed that objects fit in one of these cells. The size of the cells and the time elapsed between consecutive instants are parameters that can be adapted to particular cases.
To store absolute positions of all objects, every d time instants, GraCT uses a data structure based on the k 2 -tree, which is called snapshot. The distance, d, between snapshots is another parameter of GraCT. Between two consecutive snapshots, the trajectory of each moving object is represented as a log, which is an array of relative movements with respect to the previous time instant.

3DGraCT
3DGraCT proposes an extension of GraCT to three dimensions, so the space is divided into cells (small cubes) of fixed length, that form a bigger cube.
Snapshots. Each d time instants, there is a snapshot S k , where k is the time instant represented by the snapshot. These snapshots are organized as k 3 -trees. A leaf of the k 3 -tree set to 1 (i.e., a 1 in the bitmap L) means that one or more aircrafts 3 are placed in the corresponding cell, but the snapshot needs to determine which objects are located in that cell. Following the order of 1s in L, an array of object identifiers (aircrafts) holds that information. This array is denoted as perm, since it is a permutation [11]. An additional bitmap, called Q, is aligned with perm. It marks with 0 that the aligned object identifier in perm is the last object in the corresponding cell, and 1 means that more objects are located in that cell. Figure 3 shows an example of snapshot. 4 The two matrices models the first two slices of an 8 × 8 × 8 cube representing the 3D space. Each slice contains the horizontal positions of all aircrafts flying at a given altitude. Each matrix shows object identifiers at certain positions, and the corresponding k 3 -tree encodes this information by assuming that no objects are contained in the remaining slices. Each non-empty position in matrix corresponds to a bit set to 1 in L. The object identifiers corresponding to the first 1 in L (which is at position 3 of L) are stored starting at position 1 of perm. Q is then accessed to count the number of objects that are located in this cell: a sequential search is performed from Q[1] until the first 0 (located at Q [2]). Thus, there are two objects in the inspected cell. The corresponding object identifiers are retrieved from perm[1]=3 and perm [2]=6. Now, in position 3 of perm starts the object identifiers corresponding to the second 1 in L, and so on.
These structures allow 3DGraCT to address two types of queries: -Find the objects in a box of the 3D space. The k 3 -tree is traversed from the root to the leaves to obtain positions n 1 , n 2 , . . . n m , in L, that corresponds to positions marked with 1 in the queried box. For each n i , we count the number of 1s in the array of leaves L until the position n i ; it obtains the number of non-empty leaves up to the n th i leaf, x = rank 1 (L, n i ). Then, the position of the (x − 1)th 0 in Q is obtained, which indicates the last bit of the previous leaf (with objects), and by adding 1, we get the first position in perm with the objects of the leaf corresponding to n i , p = select 0 (Q, x − 1) + 1. From p, object identifiers aligned with 1s in Q are retrieved, until a 0 is reached (it marks the last object identifier located in a leaf).
-Find the position in the 3D space of a given object. The desired object identifier is first searched in perm. Our permutation is enhanced with shortcuts to avoid sequential searches. Assuming the object identifier is located at position k, the following step looks for its corresponding position in L. We calculate the number of leaves before the object at perm[k]: y = rank 0 (Q, k − 1). Then we find in L the position of the (y + 1) th 1, that is, select 1 (L, y + 1). This value is used to traverse the k 3 -tree upwards in order to obtain the cell position in the 3D space, and thus the horizontal position and altitude of the object.
Log of relative movements. The use of a snapshot for encoding each time instant would consume too much space, instead, between snapshots, 3DGraCT stores for each aircraft the relative movements with respect to the last known position. A relative movement consists of 3 values, x, y, z , which are the number of cells of difference between the new position and the last known position, in each dimension. Probably, x, y, z will be numbers with a small magnitude, as the differences between consecutive time instants cannot be very big. Instead of using 32 bits for each value, we fit the three values into a 32-bit integer using 12 bits for the x and y values and 8 bits for the z component. In Figure 4(a), we can see a relative movement of 1 cell up on the y-coordinate, 3 cells to the right on the x-coordinate and 2 cells down on the z-coordinate. Below, observe that those values are encoding using Zig-Zag encoding (−1 → 1, 1 → 2, −2 → 3, . . .), and then they are packed in a 32-bit integer.
Obviously, this works well as long as the assumption that there are small differences between two consecutive positions is maintained. However, there may be periods of time without information about the positions of the aircraft (for example, the aircraft is in an area without reception stations). In those cases, the 32-bit integer comprising x, y, z would not be enough. Observe that, to save space, our method does not explicitly store the time instant of a recorded position, it can be derived from its position inside the log. Therefore, 3DGraCT requires a method to manage that disappearances/appearances.
Between two consecutive snapshots S k and S k+d , each object is represented by a log, L k,k+d (id j ), where id j is the identifier of the object. It is a sequence of codewords of the following types: (1) an integer encoding a relative movement; (2) Disappearance (D) codeword, which means that we have no information about the position object id j from one time instant of L k,k+d (id j ) until its end; (3) Absolute appearance (AA), which means that we have no information about the position of id j from the beginning of L k,k+d (id j ) until a time instant covered by L k,k+d (id j ), where that information appears; (4) Relative disappearance, which means that the information about the position of id j disappears in a time instant of L k,k+d (id j ), but reappears in a time instant of the same portion of the log.
In order to maintain the synchronization of the sequences of values in L k,k+d (id j ), the appearances and disappearances require the storage of their corresponding time instant. In addition, they also require the storage of the absolute position of the appearance/disappearance. The relative disappearances imply the storage of the number of time instants they lasted and the relative movement with respect to the last known position.
In Figure 4(b), it is shown an example. The relative movements are depicted with the three relative displacements z, x, y 5 . The array D stores the duration of a relative disappearance and the exact time instant of absolute appearances and disappearances. For example, in L 0,4 (1), there is a relative disappearance that lasts two instants, and in L 0,4 (7), the object appears at time instant 3. In addition, array P stores the relative movements of relative disappearances and the absolute position of absolute appearances or disappearances. For example, in L 0,4 (1), the 1, 4, 1 tuple in P 0,4 (1) means that the object reappeared 1 cell upwards in the z-coordinate, 4 cells to the right in the x-coordinate, and 1 cell upwards in the y-coordinate. In L 0,4 (7), the object appears in the absolute position (0, 5, 2) (see P 0,4 (7)). In the figure, the values are aligned to their corresponding time instants, but this is only for illustration purposes, thanks to the array D, for one object, all the logs are stored as a sequence. D and P are compressed with DACs [3], a compressor for sequences of integers that provides direct access to any position without the need of decompressing the previous numbers.
Compressing the log. Logs represent an important saving in space with respect to snapshots, but it is possible to obtain additional compression taking advantage the following fact: aircrafts spend most of the time following the same course at a constant speed. This situation will be represented in the logs as sequences of repetitive numbers, that is, the same relative displacements with respect to the previous time instant. These series of similar numbers are compressed very efficiently using a grammar compressor, such as Re-Pair.
To improve the query processing, the Re-Pair rules in 3DGraCT are enriched with additional information. Each rule in R has the following information: s → a, b, #t, x, y, z, M BB, where: (1) s, a and b are the components of a normal rule of Re-Pair, (2) #t is the number of instants covered by the rule, (3) z, x, y are the relative coordinates of the final position of the object after the application of the rule (that is, the displacement considering (0,0,0) the initial position before the application of the rule) and, (4) M BB is the Minimum Bounding Box enclosing the movements of the rule. MBB is represented by six coordinates (z1, x1, y1, z2, x2, y2), which are the points at the ends of a diagonal of the box.
For example, in Figure 4, in L 0,4 (4), the two 1, 1, 1 consecutive relative movements produce a rule, W → 1, 1, 1 , 1, 1, 1 , 2, 2, 2, 2 , (0, 0, 0, 2, 2, 2), and then L 0,4 (4) = W, D. Thanks to the additional information, the non-terminal symbols of the logs do not need to be decompressed in many cases. For example, if we wish to know the position of object 4 at t 2 , we obtain its absolute position in the snapshot S 0 (Figure 3), which is (0,4,0), and then the first symbol of L 0,4 (4) (W ) is applied. Since W covers 2 time instants, its application to the position at t 0 produces the position of the object at the queried time instant. For this, the relative displacement (2, 2, 2) is added to the original position, obtaining the position (2, 6, 2).
Obtain the position of an object. To obtain the position of an object at a given time instant t q , first, the algorithm retrieves the position of the object in the closest snapshot to t q . If the snapshot does not represent t q , then the algorithm follows the movements through the log until it reaches t q , as it was explained in the previous section, using the relative coordinates included in the rules when possible. When the nearest snapshot is located before t q , the process follows a forward traversal of the log, otherwise, the process performs a backward traversal.
Obtain the trajectory of an object. Given an interval of time [t s , t e ] and an object, this query obtains all the positions of the object between t s and t e . First, the query obtains the position of the object at t s using the algorithm explained for the previous query, and then it applies the movements of the log until it reaches the position at t e . Since the additional information of the rules does not contain the detailed positions of the trajectory, the algorithm has to decompress every non-terminal value of the log containing a t i ∈ [t s , t e ].
Time slice query. Let r = [x 1 , y 1 , z 1 ] × [x 2 , y 2 , z 2 ] be a rectangular cuboid (or box) and t q a time instant, this query returns all objects within r at t q . Let (s x , s y , s z ) be the maximum speed vector of any object in our dataset, that is, the maximum speed in each of the three axes of the space achieved by any object in the dataset. We denote E r (t k , t q ), the expanded region of r from t k to t q , as the area that contains any object active at t k capable of being located within r at t q . Hence, E r (t k , t q ) is r extended in the three dimensions; in the x-axis to the coordinates [x 1 − s x · (t q − t k ), x 2 + s x · (t q − t k )], and repeat the same for the y-axis and z-axis. Assuming that the closest snapshot is S k and that t k ≤ t q , the algorithm obtains the candidate objects C inside E r (t k , t q ) at t k . If t k = t q , the algorithm returns C. Otherwise, it tracks the movements in L k,k+d for each object in C until it reaches t q . During this process, after obtaining the position of an object c j at t i , we can discard c j if it is outside E r (t i , t q ). The position at t q can be given by a terminal or a non-terminal value. In the first case, we apply the movement and check if the object is within r. In the second case, the object is part of the solution when the MBB of the additional information of the rule defining the non-terminal value is completely contained in r, and the object can be pruned if its MBB does not intersect r. However, when the MBB intersects r (but it is not completely contained), the algorithm has to decompress the nonterminal symbol using the Re-Pair rule to obtain the exact position of the object at t q . If the closest snapshot to t q is after it, then the algorithm performs the same process backwards.
Time interval query. Given a box r and an interval of time [t s , t e ], this query obtains all objects within r at any t i ∈ [t s , t e ]. This query could be solved in a similar way to the previous one. However, to avoid large expanded regions, that lead to track too many candidate objects, the query interval [t s , t e ] is divided into as many queries as portions of log overlaps. Then, each one of these portions [t s , t e ] can be solved in a similar way to time-slice. First, the algorithm obtains the candidates from the closest snapshot, using the expanded region with respect to t e ; then it applies the movements of the log. During the processing of the log of a candidate object c j , the algorithm has to take into account that when the traversal reaches a symbol s m that after its application obtains the position t e ), then c j can be discarded of the processing of the current portion; (3) if (x i , y i , z i ) is outside r but within E r (t i , t e ), then c j continues as a candidate that needs to be tracked. If s m is a non-terminal symbol that produces a position at t i > t e and covers the time interval [t u , t i ], where t u ≤ t e : (1) if the MBB of s m is fully within r, then c j is part of the solution (2) if the MBB of s m does not intersect r, then c j is discarded in the processing of the current portion. (3) if the MBB of s m intersects r, the algorithm has to decompress s m to check if s m involves any t l ∈ [t u , t e ] whose position is within r.

Experimental Evaluation
Our experiments analyze space/time tradeoffs of 3DGraCT using real-world ADS-B data. We also evaluate the use of interpolation to fill in large periods of missing data during the trajectory. For comparison purposes, we propose a baseline including the MVR-tree [17], but we do not include ADS-BI [21] because it does not provide altitude-based queries, and its inner index stores some string dimensions which are not covered by 3DGraCT.
Dataset details. We use four real ADS-B datasets including descriptive data of flights between different airports of Europe (see details in Appendix A). Each dataset covers a different period of time, namely one day, one week, two weeks, and one month. Positions are discretized into a cube where the cell size is 5 kilometers in x-axis, 5 kilometers in y-axis, and 100 meters in z-axis. Since aircraft positions can contain incorrect information and they can be emitted at different time rates, we discard incorrect positions and normalize timestamps to obtain regular instants every 15 seconds.
Gate-to-gate trajectories are difficult to reconstruct from ADS-B data because some broadcasted positions are lost, mainly due to lack of coverage. Although we use disappearance and reappearance codewords to represent these situations, we consider relevant to understand how they affect to 3DGraCT tradeoffs. We use the original datasets to generate a new ones, where aircraft positions are interpolated when no information is available during, at least, 15 minutes. As consequence, we have eight datasets: four real-world datasets (1D, 1W, 2W, 1M ) and four interpolated datasets (1D-I, 1W-I, 2W-I, 1M-I ). Table   Dataset   1D  1D-I  1W  1W-I  2W  2W-I  1M  1M-I  Time  1  1 shows the details of each dataset. Note that the fourth and fifth rows give, respectively, dataset sizes of binary and p7zip-compressed representations.
Compression ratio. We define compression ratio as the ratio between the binary size and the compressed size. The last row of Table 1 gives compression ratios reported by p7zip for all datasets, while Figure 5(a) illustrates 3DGraCT numbers for one day and one month datasets, using different periods of snapshot (120, 240, 360 and 720 time instants). p7zip report stable ratios around 22-24%, but 3DGraCT effectiveness is clearly influenced by the distance between snapshots, because snapshot encoding requires more space than log compression. Thus, the more-distanced the snapshots are, the better the results are. In our experiments, 3DGraCT reports its best ratios using a separation of 720 time instants between snapshots, outperforming p7zip in all datasets. For instance, 3DGraCT reports 22.29% for 1D and p7zip 23.41%. This gap increases for larger datasets: 3DGraCT only needs 14.73% of the original 1M size, while p7zip demands 22.01%. Thus, 3DGraCT is more effective than a powerful compressor like p7zip, while retaining search capabilities. This comparison also applies for interpolated datasets. Note that, in this case, 3DGraCT reports slightly better results, meaning that missing information adds an small overhead (≈ 2%) to our structure.  Figure 5(b) shows that query times of Object t increase with distance between snapshots because larger log portions must be processed. On the contrary, Figure 5(c) shows that Trajectory queries are slowler for less distanced snapshots because more snapshots must be checked.
In region queries, Time Slice and Time Interval, the number of candidates depends on the period between snapshots. Time Slice is slower as the distance between snapshot gets larger (see

5(f) and 5(g) show that
Time Interval queries behave similar to Time Slice ones, except in the right part of Figure 5(g). In this case, the expanded region covers the whole space for each period of snapshot, so the number of candidates between different settings remains constant. Thus, traversing the log demands the same computation, but less snapshots are checked for larger periods. Finally, it is worth noting that the effect of interpolation is not very relevant to 3DGraCT performance. It is only a slight improvement for region queries and large datasets. Thus, we conclude that the interpolation of missing positions avoids the cost of managing appearances and reappearances, improves Re-Pair effectiveness, and allows logs to be processed faster. For this reason, querying real-world datasets are 3%-10% slower.
Comparison with MVR-tree. 3DGraCT and MVR-tree are compared over the real-world datasets of our setup: 1D, 1W, 2W and 1M. It is worth noting that MVR-Tree space requirements are 250-300 times larger than 3DGraCT one, but we tune MVR-Tree to run on main-memory.
Our analysis show that MVR-tree is only efficient for Time Slice, Time Interval, and knn queries. Although MVR-tree can obtain the position of an object at a given time instant, or can follow the trajectory of the object in a given interval, these are expensive queries.
MVR-tree can be enhanced with an auxiliary 3DR-tree [17], but the resulting structure would consume even more space. Thus, we only analyze queries where MVR-tree is efficient. Figures 5(d) and 5(e) show that MVR-tree outperforms 3DGraCT in Time Slice queries. However, our structure is better in Time Interval queries for large intervals ( Figure 5(g)). We study the turning point where the 3DGraCT starts to improve the MVR-tree, by increasing the time interval length. Figure 5(h) shows this comparison for the 1M dataset, and a period of snapshot of 720. 3DGraCT outperforms MVR-tree for time intervals over 550 and 200 time instants in small and large regions, respectively.

Conclusions
This paper introduces 3DGraCT, a new data structure capable of compressing and querying 3D trajectories with no prior decompression. 3DGraCT extends an existing 2D compact data structure (GraCT) to support a third dimension, enabling object positions to be enhanced with descriptive altitude data. Our improvements to GraCT are more than just improving object descriptions because 3DGraCT also enables for resolving altitude-based queries.
3DGraCT has been evaluated using real-world trajectories reconstructed from ADS-B descriptions. 3DGraCT reports better compression ratios than universal compressors like p7zip (3DGraCT uses up to 50% less space), while retaining search capabilities. Compared to traditional spatio-temporal solutions, 3DGraCT needs 2 orders of magnitude less space than MVR-tree, being competitive in query performance. Finally, we also study the effect of missing subtrajectories, concluding that interpolation is effective in different cases.


Dynamic Chain Graph Models for Ordinal Time Series Data

This paper introduces sparse dynamic chain graph models for network inference in high dimensional non-Gaussian time series data. The proposed method parametrized by a precision matrix that encodes the intra time-slice conditional independence among variables at a fixed time point, and an autoregressive coefficient that contains dynamic conditional independences interactions among time series components across consecutive time steps. The proposed model is a Gaussian copula vector autoregressive model, which is used to model sparse interactions in a high-dimensional setting. Estimation is achieved via a penalized EM algorithm. In this paper, we use an efficient coordinate descent algorithm to optimize the penalized log-likelihood with the smoothly clipped absolute deviation penalty. We demonstrate our approach on simulated and genomic datasets. The method is implemented in an R package tsnetwork.

Introduction
Graphical models are an efficient tool for modeling and inference in high dimensional settings. Directed acyclic graph (DAG) models, known as Bayesian networks (Lauritzen, 1996), are often used to model asymmetric cause-effect relationships. Models represented by undirected graphs are used to model symmetric relationships, for instance gene regulatory networks.
Some graphical models are able to represent both asymmetric and symmetric relationships simultaneously. One such model so-called chain graph model (Lauritzen, 1996, Lauritzen andWermuth, 1989) which is a generalization of directed and undirected graphical models. Chain graph models contain a mixed set of directed and undirected edges. The vertex set of a chain graph can be partitioned into chain components where edges within a chain component are undirected whereas the edges between two chain components are directed and point in the same direction. Recently, chain graph models are considered in a time series setting (Abegaz and Wit, 2013, Gao and Tian, 2010, Dahlhaus and Eichler, 2003. There is a rich literature on reconstructing undirected graph for continuous data, categorical data, and mixed categorical and continuous data (Behrouzi and Wit, 2017, Mohammadi et al., 2015, Dobra et al., 2011, Hoff, 2007 and similarly for directed acyclic graphs (Colombo et al., 2012, Kalisch andBühlmann, 2007). Recently, Abegaz and Wit (2013) have proposed a method based on chain graph model for analyzing time course continuous data, like gene expression data. However, many real-world time series data are not continuous, but are categorical or mixed categorical and continuous. Until now constructing dynamic networks for non-continuous time series data has remained unexplored. Here, we develop a method to explore dynamic or delayed interactions and contemporaneous interactions for time series of categorical data and time series of mixed categorical and continuous data.
The proposed method is based on chain graph models, where the ordered time steps build a DAG of blocks and each block contains an undirected network of variables under consideration at that time point. The method developed in this paper is designed to analyze the nature of interactions present in repeated multivariate time series mixed categorical and continuous data, where we use time series chain graphical models to study the conditional independence relationships among variables at a fixed time point and "causal" relationship among time series components across consecutive time steps. The concept of causality that we use is the concept of Granger causality (Granger, 1969), which exploits the natural time ordering to achieve a "causal" ordering of the variables in multivariate time series. The idea of this causality concept is based on predictability, where one time series is said to be Granger causal for another series if the latter series to be better predicted using all available information than if the information apart from the former series had been used. Our inference procedure not only enforces sparsity on interactions within each time step, but it also between time steps; this feature is particularly realistic in a real-world dynamic networks setting.
We proceed as follow: in section 2, we explain the method where we first introduce dynamic chain graph models in section 2.1, then we propose the Gaussian copula for mixed scale time series data in section 2.2. In sections 2.3 and 2.4 we define a model for underlying multivariate time series components and we explain the procedure of penalized inference based on the L1 norm and smoothly clipped absolute deviation (SCAD) penalty terms. In section 2.5 we present a method for obtaining the loglikelihood of the observed mixed scale time series component under the penalized EM algorithm and we proceed with model selection for tuning the penalty terms. In section 3 we study the performance of the proposed dynamic chain graph model under different scenarios. Furthermore, we compare its performance with the other available methods. The proposed method is demonstrated in section 4 to investigate the course of depression and anxiety disorders.

Dynamic chain graph models
A chain graph is defined as G = (V, E) where V is a set of vertices (nodes) and E is a set of ordered and unordered pairs of nodes, called edges, which contains the directed and undirected interactions between pairs of nodes. A dynamic chain graph model is associated with a time series chain graph model, where the dependence structure of the time series components can be divided into two sets: intra time-slice dependencies, which are represented by undirected edges that specify the association among variables in a fixed time step, and a set of inter time-slice dependencies, which are represented by associations among variables across consecutive time steps. Links across time steps are directed pointing from a set of nodes at a previous time step, V (t−1) , to nodes at the current time step, V t . The dynamic chain graph models in our modeling framework relates the time series components at time t to only that of at time t − 1, but this can be easily extended to a higher order (d ≥ 2) time steps.
Let Y(t) = (Y 1 (t), . . . , Y p (t)), t = 1, . . . , T be an p-dimensional time series vector representation of p variables that have been studied longitudinally across T time points. Each time series component Y (t) is assumed to be sampled n times. Thus, Y ij (t) represents the value of the j-th variable at time t for the i-th sample, i = 1, . . . , n, j = 1, . . . , p.
Here, we focus on non-Gaussian multivariate time series data such as ordinalvalued time series taking values in {0, 1, . . . , (c k − 1)}, where c k is the number of possible categories, or mixed categorical and continuous time series data, which routinely occurs in real world settings.

Gaussian Copula
To model dependencies among p-dimensional vector y we use a Gaussian copula, defined as where Φ p (.|Ω) is the a p-dimensional Gaussian cdf with correlation matrix Ω p×p , and y = (y 1 , . . . , y p ). From equation (1) the following properties are clear: the joint marginal distribution of any subset of Y has a Gaussian copula with a correlation matrix Ω and univariate marginals F j . The Gaussian copula can be expressed in terms of a latent Gaussian variable Z = Z 1 , . . . , Z p as follow Since the marginal distributions F j are nondecreasing, observing y i 1 j < y i 2 j implies z i 1 j < z i 2 j . This can be written as set A(y) where given the observed data y j = (y 1,j , . . . , y n,j ), the latent samples z j = (z 1,j , . . . , z n,j ), are constrained to belong to the set A(y) = {z ∈ R n×p : max{z s,j : y s,j < y r,j } < z r,j < min{z s,j : y r,j < y s,j }} If an observed value of y j is missing, we define the lower bound and the upper bound of z (r) j as −∞ and ∞, respectively.

Model definition
We assume a stable dynamic chain graph model meaning that the structure of interactions within each time point remains stable for previous and current time step, and interactions between consecutive time steps are stable too. We use a vector autoregressive process of order 1, VAR(1), to describe the directed latent interactions, where t ∼ N (0, Θ −1 ) describes the undirected instantaneous interactions. The parameter set of this model contains all the conditional independence relationships in the dynamic chain graph model where the following terms hold: θ jj = 0 if and only if Z Given the set A(y), we calculate the likelihood as Given the set of parameters, the event z ∈ A(y) in (4) does not depends on marginals and contains the relevant information about the copula and the parameters of interest Θ and Γ.
We drop the second term in (4) because this term does not provide any information about intra and inter time-slice dependencies. As Hoff (2007) proposes we use f Z (z ∈ A(y) | Θ, Γ) as the rank likelihood, We ignore the second term in (5) as we do not want to make additional assumption on the unconditional distribution of Y (1) . And we start from t = 2, where we compute the conditional log-likelihood using the conditional distribution f (z (t) |z (t−1) ). According to (3) the conditional distribution Z (t) | Z (t−1) follows a multivariate normal distribution which its density for t-th observation is defined as

Penalized EM inference
In Gaussian copula, we treat the marginals distributions as nuisance parameters since our main goal is to learn the dependence structure among time series components both at a fixed time step t ∈ N and also across consecutive time steps. We use an empirical 1 n 1(y ij ≤ y) (Genest et al., 1995) to estimate marginals. Genetic time series data often are high dimensional due to a large number of variables that are measured on small number of samples across only few time steps. Furthermore, many real-world networks (e.g. genetic, genomics, and brain networks) are intrinsically sparse. Thus, incorporating sparsity into the proposed dynamic chain graph model makes the derived model more biologically plausible. Accordingly, we propose a dynamic chain graph model for genetic data based on the penalized likelihood. In order to find the penalized maximum likelihood estimation we will use the EM algorithm (Green, 1990). This modeling technique provides sparse estimates of the autoregressive coefficient matrix Γ and the precision matrix Θ in (3) which are used to reconstruct inter and intra time-slice conditional independences, respectively.
The E-step of the EM algorithm is given by Under the assumption described in (6), the E-step can be written as such that conditional expectation at current time, S cc , and at past, S pp , is defined as and the conditional expectation at inter time-slice dependence is The latent variables Z i,p } is used to calculate the conditional expectation of intra time-slice dependencies S pp and S cc , respectively. And Z i,p } is used to calculate S pc . All the three above mentioned conditional expectations are a p×p matrix. When j = j they can be computed through the second moment E(Z (t) 2 ij | y i ; Θ , Γ ). When j = j we use a mean field theory approach (Chandler, 1987) to approximate them as for intra time-slice dependencies, and for inter time-slice dependencies follows as This approximation performs well when the interaction between Z (t) i,j and Z (t) i,j given the rest of the variables, and the interaction between Z (t−1) i,j and Z (t) i,j given the rest of the variables are close to be independent; this often holds in our proposed dynamic chain graph model which Θ and Γ are sparse.
When j = j the off-diagonal elements of S cc , S pp , and S pc matrices can be computed through the first moment as and the second moments is Given the property of Gaussian distribution, (Z ; Θ, Γ follows a multivariate normal distribution with mean and variance-covariance matrix (13) and (14) follows a multivariate normal distribution with mean µ ij and variance v ij as follow Calculating the exact value of the first and second moments is computationally expensive. Moreover, we approximate the first and the second moments as follow The conditional distribution of Z (15) and (16) follows a multivariate normal distribution with mean and variance-covariance matrix as follow ij +1 ] which the first and second moments can be obtained via lemma 2.1.
Both means µ i,j and µ i,j are a linear function of z i,−j . Applying Lemma 2.1 on the conditional expectations in (15) and (16) leads to following approximations where δ Here, the first order delta method is used to approximate the nonlinear terms [more details in (Guo et al., 2015)]. Moreover, we approximate the elements of inter time-slice conditional expectation matrix S pc through equations (17) and (18). For approximating the elements of intra time-slice conditional expectation matrices S pp , S cc we refer to the Appendix.
The M-step of the EM algorithm contains two-stage optimization process where we maximize expectation of the penalized log-likelihood with respect to Θ and Γ. We introduce two different penalty functions P λ (.) and P ρ (.) for intra time-slice conditional independencies Θ, and inter time-slice conditional independencies Γ, respectively. Therefore, the objective function for optimization can be defined as denotes the expectation of S Γ given the data and updated parameters, and θ jj and γ jj are the jj -th element of the Θ and Γ matrices. Among different penalty functions, we consider the L 1 norm and smoothly clipped absolute deviation (SCAD) penalty functions which have the desirable sparsity properties.
The L 1 penalty leads to a desirable optimization problem, where the log-likelihood is convex and can efficiently be solved using various optimization algorithms at the k-th iteration of the EM. Under this penalty function, the updated estimates are given via where the sparsity level of intra and inter time-slice conditional independences are controlled by λ and ρ. L 1 penalty is biased due to its constant rate of penalty.
To address this issue, Fan and Li (2001) proposed SCAD penalty, which results in unbiased estimates for large coefficients.
SCAD penalized EM. The SCAD penalty function is expressed as where λ and a are two tunning parameters. The function P λ,a (θ) corresponds to a quadratic spline on [0, ∞) with knots at λ and aλ. A similar function can be written for P ρ,a (γ) where ρ and aρ are two knots. The SCAD penalty is symmetric but non-convex, whose first order derivative is given by P λ,a (θ) = λ I(|θ| ≤ λ) + (aλ − |θ|) + (a − 1)λ I(|θ| > λ) , a > 2 The notation z + stands for the positive part of z. Fan and Li (2001) showed that in practice a = 3.7 is a good choice. Maximizing non-convex penalized likelihood is challenging. To address this issue, we use an efficient algorithm proposed in Fan et al. (2009), which is based on local linear approximation, to maximize the penalized loglikelihood for the SCAD penalty function. In each its step, a symmetric linear function is used to locally approximates the SCAD penalty. Using the Taylor expansion, P λ,a (θ) and P ρ,a (γ) can be approximated in the neighbor of θ 0 and γ 0 as follow: Due to the monotonicity of P λ (.) and P ρ (.) over [0, ∞), the derivatives P λ (.) = ∂ ∂θ (P λ (θ)) and P ρ (.) = ∂ ∂γ (P ρ (γ)) are non-negative for θ ∈ [0, ∞) and γ ∈ [0, ∞). Therefore, under the penalized log-likelihood with SCAD penalty, the estimate of the sparse parameters Θ (k) and Γ (k) relies on the solution of the following optimization problem at step k jl are jj -th element of Θ and jlth element of Γ, respectively. The SCAD penalty applies a constant penalty to large coefficients, whereas the L 1 penalty increases linearly as |θ| increases. This features keep the SCAD penalty against producing biases for estimating large coefficients. Therefore, the SCAD penalty overcome the bias issue of the L 1 penalty. Then a two stage-optimization problem within the M-step of the EM algorithm is employed to solve the objective functions (20) or (21) to estimate the parameters Θ and Γ.
Glasso calculation of Θ (k) . For the SCAD penalty-based estimation, in the first stage we optimize for previous Γ . This optimization can be solved efficiently using the graphical lasso algorithm proposed by Friedman et al. (2008). Due to the sparsity in each iteration, we consider a one-step local linear approximation algorithm (LLA). Zou and Li (2008) showed that one-step LLA, asymptotically, performs as well as the fully iterative LLA algorithm as long as initial solution is good enough. In practice, we take the initial value as the L 1 penalty graphical LASSO for estimating the intra time-slice conditional independences Θ in order to calculate the initial weights w jj and ν jl .
Regularized coordinate descent algorithm for Γ (k) . After we finish an updating Θ in the first-stage of the optimization, in the second-stage we proceed to update the estimate of Γ given the updated Θ. In the SCAD penalty-based we optimize This objective function is quadratic in Γ for given Θ (k) λ . Thus, we use a direct coordinate descent algorithm to calculate Γ (k) ρ . So, the derivative of the penalized negative log-likelihood (22) with respect to γ jl is where sgn(.) is the sign function. These are the Karush-Kuhn-Tucker (KKT) equations defining the solution to the maximization problem. We note that for an arbitrary matrix A p×p , ∂tr(ΓA)/∂γ jl = a lj = e l Ae j , where e l and e j are the corresponding base vector with p dimension each. Setting the derivative of negative log-likelihood (23) to zero, we get an update for the elements of Γ matrix as follow λ )e j }, γ jl , and Γ (k) ρ are the estimates in the last step of the iteration inside the optimization (24).
Given the two-stage optimization problem inside the M-step, we update the S Γ matrix in the E-step. This iterative procedure continues until the difference between previous (Θ ρ ) becomes smaller than a, user specified, tolerance. Based on our simulation experiments, the EM algorithm converges in a few iterations (at most 5 iterations is needed to reach the convergence). We define the estimate as the stationary point of the EM, ( Θ λ , Γ ρ ) = lim Table 1 Performance measure results of the simulation study for tsnetwork and SparseTSCGM using SCAD penalized likelihood estimation for the precision and autoregressive coefficient matrices for fixed time point, t=5. In SparseTSCGM* the normal transformation is applied to the simulated ordinal data.

Selection of tuning parameters
To determine the sparsity of the proposed dynamic chain graph model, the tunning parameters λ and ρ have to be tuned. We focus on estimating the sparse intra and inter time-slice conditional independences Θ and Γ, we employ the Bayesian information criteria (BIC) to select the tuning parameters λ and ρ, where T and p are the number of time points and the number of variables, respectively, and df( Θ λ ) shows the number of non-zero elements in the off-diagonal of Θ λ , and df( Γ ρ ) is the number of non-zero elements of Γ ρ . The approximation made in BIC is the result of a Laplace-type of approximation, which makes fast calculation feasible. We choose the optimal value of the penalty parameters that minimizes BIC(λ, ρ) on a grid of candidate values for λ and ρ. One may consider other information criteria that suits for graph estimations. Wang et al. (2007) and Yin and Li (2011) has been shown that BIC performs well for selecting the tunning parameter of penalized likelihood estimation.

Simulation study
To investigate and assess the performance of the proposed dynamic chain graph model, we set up a simulation to generate sparse Θ and Γ matrices similar to Abegaz and Wit (2013), and Yin and Li (2011). Here we evaluate the performance of the proposed method with respect to different random graph structures for Θ and Γ matrices. Simulating different graph structures for Θ can be performed through the R package flare. For generating Γ matrix we took the upper diagonal of an independently generated Θ a long with a 0.2% nonzero diagonal elements sampled from uniform (0, 1), similar to the R package SparseTSCGM. Table 2 Performance measure results of the simulation study for tsnetwork and SparseTSCGM using SCAD penalized likelihood estimation for the precision and autoregressive coefficient matrices for fixed time point, t=10. In SparseTSCGM* the normal transformation is applied to the simulated ordinal data. First we simulate data from N p (0, Θ −1 ) at time t = 1, for the next time steps t = 2, . . . , T we use VAR(1) model such that Z (t) |Z (t−1) ∼ N (ΓZ (t−1) , Θ −1 ). Then, n i.i.d samples is generated for each time point. This results in p-variate time series data. Finally, we discretize the obtained time series data with Gaussian marginals into randomized quantile ranges and treat them as categorical time series data. The simulations are repeated 50 times independently for different values of p, n, t.
To assess the performance of our proposed method in recovering the intra and inter conditional independence relationships we compute the F 1 -score, sensitivity and specificity measures, which are defined as: where TP, TN, FP, and FN are the numbers of true positive, true negative, false positive, false negative in identifying the non-zero elements in the Θ and Γ matrices. There are four categories in NESDA data: (i) sleep in blue, (ii) mood in green, (iii) appetite in yellow, (iv) somatic in gray, (v) mental in red.
We note that high values of the F 1 -score, sensitivity and specificity indicate good performance of a method for the given combination of p, n and t. However, as there is a natural trade off between sensitivity and specificity, we focus particularly on the F 1 -score to evaluate the performance of each method. We compare the finite sample performance of the proposed approach using SCAD penalized maximum likelihood with a recently proposed approach implemented in R package SparseTSCGM (Abegaz et al., 2015). For further comparison we have applied SparseTSCGM to the original simulated ordinal data and to the transfered data using the normal transformation. We present the simulation results of sparse precision and autoregressive coefficient matrices in Table 1 and Table 2 based on optimal tuning parameters chosen by the minimum EBICs. In each simulation setting, we have very sparse matrices with only (1/p) × 100 nonzero entries. From the tables, we can see that in most cases our method scores better in terms of the F1-score compare with the alternative method. These results suggest that, though recovering sparse network structure in ordinal time series data is a challenging task, the proposed approach has a good performance on model-based simulations. We note here that improved model performance can be gained by allowing the tuning parameters ρ and λ to vary with each simulation.

Netherlands Study of Depression and Anxiety
We applied our method to a Netherlands Study of Depression and Anxiety (NESDA) Severity of Depression dataset. Depression and anxiety disorders are common at all ages. Approximately one out of three people in the Netherlands will be faced with one of these disorders at some time during their lives. It is still not clear why some people recover quickly and why others suffer for long periods of time. The Netherlands Study of Depression and Anxiety (NESDA) was therefore designed to investigate the course of depression and anxiety disorders over a period of several years. The main aim of NESDA is to determine the (psychological, social, biological and genetic) factors that influence the development and the long-term prognosis of anxiety and depression. The data consist of the 28 items (variables) that have been collected in 3 time intervals. For each of 28 variables there are four corresponding answers 0=None, 1=Mild, 2=Moderate, 3=Severe. For example, for the item "Feeling sad" there are four corresponding answers from "0" that is indicative of no depression (e.g., "I do not feel sad") to "3" referring to a more severe depressive symptom (e.g., "I feel sad nearly all the time"). A total score is derived (possible range: 0-84), and higher scores are indicative of relatively severe depressive symptomatology. From the 1799 participants, we have selected 200 patients that have been more informative. The BIC criterion selects the penalty values λ = 0.19 and ρ = 0.23. The resulting instantaneous and delayed interaction network among the 28 items are shown in Figure 1, left and right panels, respectively. Figure 1(a) shows the undirected links that suggest contemporaneous interactions among 12 items and Figure 1(b) displays the directed edges that indicate grangercausality relationships or delayed interactions between these 12 items. It is observed that item "Feeling sad" is the hub in both figures, suggesting that it plays a fundamental role in treating depression and anxiety disorders. Also, Figure 1(b), shows that there are several directed links pointing from mood category to mental category suggesting that mood disorders influence the development of mental disorders in long term. Interestingly, Figure 1b shows that sleeping disorders do not have any effect on other symptoms of depression.

Discussion
We have presented a dynamic model for multivariate ordinal time series data which assumes a chain graph representation of the conditional independence structure among time series components. The proposed model combines the Gaussian copula graphical models and dynamic Bayesian networks to infer instantaneous conditional dependence relationships among time series components and dynamic or delayed interactions possibly potentially "causal" relationships among variables at consecutive time steps. The directed edges reflect Granger causality whereas the contemporaneous dependence structure is represented by undirected edges.
To obtain sparse estimates for the instantaneous conditional dependence graph and for the Granger-causality graph, we considered penalized log-likelihood estimation using the L 1 and SCAD penalties. Simulation studies show that the proposed sparse estimates reflect the underlying intra-and inter-time slice conditional dependence networks more accurately compared to the only sparse alternative method.
The method was applied to the Netherlands study of depression and anxiety categorical time series data. The model does, however, have much wider applicability to any multivariate mixed continuous and discrete time series data.

Appendix
Another approximation that can be replaced in (15) and (16) follows as where Z (t) i,−j represents a set that contains all the variables at time step t except the j-th variable.