Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
62,260 characters · 17 sections · 34 citation commands
Score-Driven Exponential Random Graphs: A New Class of Time-Varying Parameter Models for Temporal Networks
\preprint{AIP/123-QED}
\email{[email removed]}
A network is a useful abstraction for a system composed of several single elements with some pairwise relations. The simplified description of social, economic, biological, and transportation systems, often very complex, in terms of nodes and links, attracted and still attracts an enormous amount of attention RevModPhys.74.47,bullmore2009complex,newman2010networks,easley2010networks,allen2009networks. Formally, a network $G$ is a pair $ \left( V, E \right) $ where $V$ is a set of nodes and $E$ is a set of node pairs named links. The nodes are labeled, and a link is identified by the pair of nodes it connects $(i,j)$. To each $G$, we can assign an adjacency matrix ${\mathbf{Y}}$ such that ${Y_{ij}} = 1 $ if link $\left(i,j\right)$ is present in $E$ and ${Y_{ij}} = 0 $ otherwise. Links may have orientations. The corresponding network is dubbed directed. If the elements of the adjacency matrix are allowed to be different from $0$ or $1$, one speaks of weighted networks. In the following, we will focus on directed networks rather than consider the weighted variant.
Often, systems that are fruitfully described as networks evolve in time. When the number of nodes and/or pairwise interactions change over time, one usually speaks of temporal networks holme2012temporal,craig2014interbank,rossetti2018community. This paper will focus on temporal networks where links evolve in discrete time. A temporal network is a sequence of networks, each associated with an adjacency matrix and observed at $T$ different points in time. The whole time series is given in terms of a sequence of matrices $\left\{{Y_{ij}^{\left(t\right)}}\right\}_{t=1}^{T}$.
We introduce an innovative approach to temporal networks based on two main ingredients: (i) a parametric probabilistic model, according to which one can sample a network realization. A natural choice is the class of statistical models for networks known as Exponential Random Graph Models (ERGMs). (ii) A simple mechanism to induce dynamics on the network sequence by introducing time variation on the model parameters. The Dynamic Conditional Score (DCS) approach provides a flexible candidate. Our extension of the ERGM framework allows model parameters to change over time in a score-driven fashion. We develop a new class of temporal network models and show its versatility and effectiveness in capturing time-varying features. The information encoded in $\mathcal{F}_{t-1}$ is exploited to filter the time-varying parameters (tvps) $\theta^{\left(t\right)}$ at time $t$. We refer to this class as Score-Driven Exponential Random Graph Models (SD-ERGMs). A generic SD-ERGM can be used either as a data-generating process (DGP) to sample synthetic sequences of graphs or as an effective filter of latent tvps, regardless of the true DGP.
We are by no means the first to discuss models for temporal networks. For a review of latent space temporal network models, one can refer to Ref. \onlinecite{kim2018review} Extensions of the ERGM framework for the description of temporal networks exist in the literature. Two are the main streams. The first one, termed TERGM, was pioneered in Ref. \onlinecite{robins2001random} and further explored in Refs. \onlinecite{hanneke2010discrete,cranmer2011inferential,krivitsky2014separable}. This approach builds on the ERGM but allows the network statistics to define the probability at time $t$ to depend on current and previous networks up to time $t-K$. This $K$-step Markov assumption is a defining feature of the TERGMs. A second approach allows for the parameters of the ERGM to be time-varying. A notable example is the Varying-Coefficient-ERGM lee2020varying, allowing for smooth parameter time variation. The approach differs from ours in several respects. Specifically, to infer parameter time-variation at time $t$, it uses all the available observations, including those from future times $t^\prime > t$ \footnote{In time series jargon, it is a smoother and not a filter.}. Consequently, it cannot be used to draw causal sequences of time-evolving networks. A related but different approach mazzarisi2020dynamic considers the possibility of a random evolution of node-specific parameters. As a crucial difference with the SD-ERGM, the parameter evolution is driven by an exogenous source of randomness. Following the language of Ref. \onlinecite{cox1981statistical}, the approach is parameter driven, while we consider an observation driven dynamics. Finally, it is important to mention that the social science literature has considered alternative frameworks for modeling temporal networks. Notable examples are the Stochastic Actor Oriented model snijders1996stochastic and the Relational Event Model butts20084. For an overview of contributions, we refer to the literature therein.
The rest of the paper is organized as follows. In Section (ref), we review some key concepts on ERGMs and observation-driven models. In Section (ref), we introduce the new class of models and validate it with extensive numerical experiments for three specific instances of the SD-ERGM. Section (ref) presents the results from an application to two real temporal networks: The e-MID interbank network for liquidity supply and demand and the U.S. Congress co-voting political network. Section (ref) draws the relevant findings.
The two main ingredients in our approach are the ERGMs, a static class of network models well-known in physics, and the DCSs, a recent development in time-series econometrics.
A statistical network model can be specified by providing the probability distribution over the set of possible adjacency matrices Kolaczyk:2009:SAN:1593430. If the distribution belongs to the exponential family barndorff2014information, then the model is named ERGM, and its log-likelihood takes the form
where $h$ are network statistics, $\theta$ is the vector of parameters whose component $\theta_s$ is associated with the network statistic $h_s\left({\mathbf{Y}}\right)$, and $\mathcal{K}\left(\theta\right) = \sum_{\left\{{\mathbf{Y}}\right\}} e^{ \theta_s h_s\left({\mathbf{Y}}\right) } .$ The ERGMs literature is vast and still growing schweinbergerexponential. The ERGM framework is intrinsically linked to the very well-known principle of maximum entropy Shannon:2001:MTC:584091.584093 and its applications to statistical physics PhysRev.106.620. Indeed, an ERGM with sufficient statistics $h\left(\theta\right)$ naturally arises when looking for the probability distribution which maximizes the entropy under a linear equality constraint on the statistics $h\left(\theta\right)$ PhysRevE.70.066117,PhysRevE.78.015101. The sufficient statistics $h_s\left({\mathbf{Y}}\right)$, known as network statistics, are functions of the adjacency matrix ${\mathbf{Y}}$, whose entries are binary random variables. The probability mass function (PMF) is defined by (ref). The normalizing factor $\mathcal{K}\left(\theta \right)$ is often unavailable as a closed-form function of the parameters $\theta $.
In the following, we will focus on two specific examples of ERGMs that describe distinct features of the network and require different approaches to parameter inference. The first one is meant to capture the heterogeneity in the number of connections each node can have, and it allows for straightforward maximum likelihood estimation (MLE) chatterjee2011random. It is known as beta model, fitness model, and configuration model zermelo1929berechnung,Holland81anexponential,PhysRevLett.89.258702,PhysRevE.70.066117,PhysRevE.78.015101,chatterjee2011random. The second one is an ERGM having as statistic the Geometrically Weighted Edgewise Shared Partners (GWESP). That is a network statistic describing transitivity in the formation of links, i.e., the tendency of connected nodes to have common neighbors and belongs to a family of network statistics referred to as curved exponential random graphs, proposed in Refs. \onlinecite{snijders2006new,robins2007recent} and discussed in Ref. \onlinecite{hunter2006inference}. In the latter case, the inference is complicated because the normalizing factor in (ref) is not available in closed form. In such cases, there are two standard approaches to ERGM inference, both consisting of maximizing alternative functions that are known to share the same optimum as the exact likelihood. In Appendix (ref), we provide more details on the beta model and GWESP ERGMs definitions as long as more details on the associated inference procedures.
The second main ingredient of this work is the class of DCS models creal2013generalized,harvey2013dynamic, also known as Generalized Autoregressive Score models \footnote{http://www.gasmodel.com/index.htm for the updated collection of papers dealing with GAS models.}. In the language of Ref. \onlinecite{cox1981statistical}, DCSs belong to observation-driven models. Let us consider a sequence of observations $\left\{y^{\left(t\right)}\right\}_{t=1}^T$, where each $y^{\left(t\right)} \in\mathbb{R}^M$, and a conditional probability density $P\left(y^{\left(t\right)}\vert f^{\left(t\right)}\right)$, that depends on a vector of tvps $f^{\left(t\right)} \in \mathbb{R}^K$. Defining the score as
a Score-Driven model assumes that the recursive relation
rules the time evolution of $f^{\left(t\right)}$, with $w$, $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ are static parameters, $w$ being a $K$ dimensional vector and $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ $K\times K$ matrices. $\boldsymbol{S^{\left(t\right)}}$ is a $K\times K$ scaling matrix usually chosen as a power of the inverse of the Fisher information matrix associated with $P\left(y^{\left(t\right)}\vert f^{\left(t\right)}\right)$.
The parameter updating rule can be intuitively motivated by general assumptions based on information theory principles. One can assume that the network behavior varies based on surprise: The more an observation of the network’s state, i.e., the adjacency matrix, is “unexpected", the more the relations between its components will change. The most common measure of surprise is minus the logarithm of the likelihood of observing the current state conditional on the level of the model parameters. As a second principle, one assumes that the reaction to surprise is to adapt to it, making what had been unexpected at that moment less surprising in the future. This implies that the parameters change to minimize the surprise, i.e., increase the log-likelihood of the last observation and thus move along the steepest direction the gradient provides. Then, the updated parameter value will be a linear combination of the current value and the log-likelihood score.
The structure of the conditional observation density determines the score, from which the dependence of $f^{\left(t+1\right)}$ on the vector of observations $y^{\left(t\right)}$ follows. When the model is viewed as a DGP, the update results in stochastic dynamics precisely thanks to the random occurrence of $y^{\left(t\right)}$. When the score-driven recursion is regarded as a filter, the update rule in (ref) is used to obtain a sequence of filtered $\left\{\hat{f}^{\left(t\right)}\right\}_{t=1}^T$. In this setting, one estimates the static parameters by maximizing the log-likelihood of the whole sequence of observations.
A second look at eq. (ref) reveals the similarity of the score-driven recursion with the iterative step from a Newton algorithm, whose objective function is precisely the log-likelihood function. As mentioned above, at each step, the score pushes the parameter vector along the log-likelihood steepest direction. Moreover, there are motivations, grounded on the variation of the Kullback-Leibler divergence, for the optimality of the score-driven updating rule blasques2015information, as we review in Appendix (ref).
Many well-known econometrics models can be expressed as Score-Driven models. Famous examples are the Generalized Autoregressive Conditional Heteroskedasticity (GARCH) model BOLLERSLEV1986307, the Exponential GARCH model nelson1991conditional, the Autoregressive Conditional Duration model engle1998autoregressive, and the Multiplicative Error Model engle2002new. The introduction of this framework in its full generality opened the way to applications in various contexts.
Before moving to the most crucial section, let us mention two relevant technical aspects for the applications discussed in the paper: the computation of the confidence bands of the filtered parameters and testing for the parameter temporal variation. We postpone the technical discussion to Appendix (ref) for brevity.
This section describes the methodological innovation introduced by the manuscript. We present the general SD-ERGM framework, discuss the score-driven approach's applicability to three different ERGMs in detail, and validate their performances with extensive numerical simulations.
We apply the score-driven methodology to ERGMs to allow any of the parameters $\theta_s$ in (ref) to have a stochastic evolution driven by the score of the static ERGM model, computed at different points in time. This approach results in a framework for describing temporal networks, more than in a single model, in the same way ERGM is considered a modeling framework for static networks.
Conceptually, applying the score-driven approach is pretty straightforward. Given the observations $\left\{{Y_{ij}^{\left(t\right)}}\right\}_{t=1}^{T}$, we can apply the update rule in (ref) to all or some elements of $\theta$, each of which is associated with a network statistic in (ref). To do this, we need to compute the derivative of the log-likelihood at every time step, i.e., for each adjacency matrix ${\mathbf{Y}}^{\left(t\right)}$. For the general ERGM, the elements of the score take the form
and the vector of tvps evolves according to (ref) with $f^{\left(t\right)}$ replaced by $\theta^{\left(t\right)}$. Hence, conditionally on the value of the parameters $\theta^{\left(t\right)}$ at time $t$ and the observed adjacency matrix ${\mathbf{Y}}^{\left(t\right)}$, the parameters at time $t+1$ are deterministic. When used as a DGP, the SD-ERGM describes stochastic dynamics because, at each time $t$, the adjacency matrix is not known in advance but must be randomly sampled from $P\left({\mathbf{Y}^{\left(t\right)}}\vert \theta^{\left(t\right)}\right)$ and used to compute the score. When the sequence of networks $\left\{{\mathbf{Y}^{\left(t\right)}}\right\}_{t=1}^T$ is observed, the static parameters $\left(w,\boldsymbol{\beta},\boldsymbol{\alpha}\right)$, that best fit the data, can be computed via MLE. Taking into account that each network ${\mathbf{Y}^{\left(t\right)}}$ is independent of all the others conditionally on the value of $\theta^{\left(t\right)}$, the log-likelihood can be written as
The computation of the normalizing factor and its derivative with respect to the parameters is essential for the SD-ERGM. Not only does it enter the definition of the update, but it is also required to optimize (ref).
Our primary motivation for introducing the SD-ERGM is to describe the time evolution of a sequence of networks using the evolution of the parameters of an ERGM. From the context or previous studies of static networks in terms of ERGM, we assume we know which statistics are more appropriate in describing a given network. Hence, we do not discuss the choice of statistics in the context of temporal networks but refer the reader to Refs. \onlinecite{goodreau2007advances,hunter2008goodness,SHORE201516} for examples of feature selection and goodness-of-fit evaluation.
In this final paragraph, we anticipate the SD extensions of ERGMs with given statistics detailed in the following pages. The first example allows for the exact computation of the likelihood, but the number of parameters can become significant for a large network. The second example discusses how an SD-ERGM can be defined when the log-likelihood is not known in closed form. Using extensive numerical simulations, we show that SD-ERGMs are very efficient at recovering the paths of tvps when the DGP is known, and the score-driven model is employed as a misspecified filter. Moreover, we show the first application of the Lagrange Multiplier (LM) test calvori2017testing in assessing the time-variation of ERGM parameters.
Our first specific example is the Score-Driven version of the beta model, introduced in Sec. (ref) and further discussed in Appendix (ref). We start with this model because of its wide applications and relevance in various streams of literature and because the likelihood of the ERGM and its score can be computed exactly. Moreover, the number of local statistics, the degrees, and parameters can become very large for large networks. Since we must describe the dynamics of many parameters, this last feature challenges any time-varying parameter version of the beta model. At the end of this Section, we will show how the SD framework allows for a parsimonious description of such high-dimensional dynamics.
As anticipated, the SD-beta model is defined by applying (ref) to each of the ${\overrightarrow{\theta}}$ and ${\overleftarrow{\theta}}$ parameters. Among the possible choices, we use as scaling the diagonal matrix $S_{ij}^{\left(t\right)} = {\delta_{ij} I_{ij}^{\left(t\right)}}^{-1/2} $, where $\boldsymbol{I^{\left(t\right)}} = {\mathbb E}[{\nabla^{\left(t\right)} {\nabla^{\left(t\right)}}^\prime}] $, i.e., we scale each element of the score by the square root of its variance. It is widespread, in score-driven models with numerous tvps, to restrict the matrices $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ of (ref) to be diagonal. In this work, we consider a version of the score update having only three static parameters $\left(w_s,\beta_s,\alpha_s\right)$ for each dynamical parameter $\theta_s$. The resulting update rule for the beta model is
where the superscripts $in$ and $out$ indicate the first and second half of the parameter vectors, respectively. To simplify the inference procedure, we consider a two-step approach. First, we fix the node-specific parameters $w_i$ to target the unconditional means of ${\overleftarrow{\theta}}$ and $ {\overrightarrow{\theta}} $ resulting from an ERGM with static parameters. Conditionally on the target values, we estimate the remaining parameters $\alpha^\text{in}$, $\alpha^\text{out}$, $\beta^\text{in}$, and $\beta^\text{out}$. We verified that the bias introduced by the two-step procedure is negligible, and results remain similar when the joint estimation is performed.
As mentioned in the Introduction, SD-ERGMs (as other observation-driven models, e.g., GARCH) can be seen as DGPs or predictive filters Nelson96 since tvps follow one-step-ahead predictable processes. In this Section, we show the ability of the ERGMs in the latter setting. Specifically, we simulate generic non-stationary evolution for temporal network parameters $\theta^{\left(t\right)}$. We then use the SD-ERGM to filter the paths of the parameters and evaluate its performances. It is important to note that the parameters' simulated dynamics differ from the score-driven ones specified for the estimation.
In practice, at each time $t$, we sample the adjacency matrix from the PMF of an ERGM with parameters\footnote{In the following, the notation with a bar refers to the true parameters used in the DGP.} $ \bar{\theta}^{\left(t\right)} $, evolving according to known temporal patterns that define different DGPs. We then use the realizations of the sampled adjacency matrices to filter the patterns. We consider a sequence of $T = 250$ time steps for a network of $10$ nodes, each with parameters $\overline{{\overleftarrow{\theta}}_i}^{\left(t\right)} $ and $\overline{{\overrightarrow{\theta}}_i}^{\left(t\right)} $ evolving with predetermined patterns. We test four different DGPs. The first one is a naive case with constant parameters $\overline{\theta}^{\left(t\right)} = \overline{\theta}_{0}$. The elements of $\overline{\theta}_{0}$ are chosen to ensure heterogeneity in the expected degrees of the nodes under the static beta model. For the remaining three DGPs, half of the parameters are static, and half are time-varying, evolving with either a deterministic sinusoidal function, a deterministic step function, or a stochastic AR(1) dynamics. More details on the definition of such DGPs are given in the Appendix (ref).
In the following, we benchmark the performance of the SD-ERGMs with that of a sequence of cross-sectional estimates of static ERGMs, i.e., one ERGM estimated for each $t$. We quantify the performance of the two approaches computing the Root Mean Square Error $ \frac{1}{T} \sqrt{\sum_t \left(\bar{\theta}_s^{\left(t\right)} - \hat{\theta}_s^{\left(t\right)} \right)^2 }$, that describes the distance between the known simulated path and the filtered. We then average the RMSE across all the tvps and $100$ simulations and report the results in Table (ref). These results confirm that the SD beta model outperforms the standard beta model in recovering the true time-varying pattern. Notably, this holds even when the DGP is inherently nonstationary, as in the case of the DGP, where each parameter has a step-like evolution. Indeed, the results of this Section and Section (ref) confirm that, while the SD update rule (ref) defines a stationary DGP creal2013generalized, using SD models as filters, we can effectively recover nonstationary parameters' dynamics.
Our last numerical simulations for the SD beta model explore its applicability and performance for networks of increasing size. We explore this setting for two reasons. The first one is that networks with a large number of nodes describe many real systems. The second reason is that we want to compare the performance of our approach with that of the standard beta model in regimes where the latter is known to perform better under suitable conditions. Indeed, as mentioned in Appendix (ref), asymptotic results on the single observation estimates chatterjee2011random guarantee that, if the network density remains constant as $N$ grows, the accuracy of the cross-sectional estimates increases. We want to check numerically that, within the regime of dense networks, the accuracy of the static and SD versions of the beta model reaches the same level. To check whether the SD approach provides any advantage for large networks, we perform numerical experiments similar to the previous ones but in a different and more realistic regime of sparse networks, i.e., keeping the average degree constant. Moreover, to ease the computational burden for the estimates, we consider a restricted version of the SD-Beta model, as detailed in Appendix (ref), having only one set of parameters $\left({\beta}^{\text{in}}, {\beta}^{\text{out}}, {\alpha}^{\text{in}}, {\alpha}^{\text{out}}\right)$ for the whole network, instead of one set per each node.
This analysis considers only one dynamical DGP and many different values of $N$. Among the DGPs used above, we focus on the one with smooth and periodic time variation. Most importantly, we set a maximum degree attainable for a node and let it depend on $N$ in two distinct ways, each corresponding to a different density regime: one generating sparse networks and the other dense ones. It is worth noticing that the asymptotic results of Ref. \onlinecite{chatterjee2011random} are expected to hold only in the dense case.
The average densities for different values of $N$ in the two regimes are shown in the left panel of Figure (ref). Then, for both regimes and each value of $N$, we compute the average RMSE across all tvps and all Monte Carlo replicas. In the right panel of Figure (ref), the average RMSEs for different values of $N$ indicate that, also for large networks, the SD version of the beta model attains better results compared with the cross-sectional estimates. As expected, both approaches reach the same accuracy in the dense network regimes as long as $N$ becomes larger. However, in the more realistic sparse regime, the performance of the SD-ERGM remains much superior for both small and large network dimensions.
As mentioned earlier, the dependence of the normalizing function on the $\theta$ parameters is often unknown. This fact prevents us from computing the score function and directly applying the update rule (ref) to a large class of ERGMs. To circumvent this obstacle, instead of the unattainable score of the exact likelihood, we propose to use the score of the pseudo-likelihood, discussed in Sec. (ref), that we refer to as pseudo-score
in place of the exact score in the definition of the SD-ERGM update (ref). Additionally, we use the pseudo-likelihood for each observation ${\mathbf{Y}}^{\left(t\right)}$ in (ref) to infer the static parameters.
Our approach, based on the score of the pseudo-likelihood, requires as input the change statistics for each function $h_s\left({\mathbf{Y}^{\left(t\right)}}\right)$ \footnote{For practical applications, it is very convenient that, for a large number of network functions, an efficient implementation to compute change statistics is made available in the R package ergm hunter2008ergm.}. In the following, we show that the update based on the pseudo-likelihood score effectively filters the path of tvps. Remarkably, this is true even when the probability distribution in the DGP is exact, i.e., when we sample from the exact likelihood and then use the SD-ERGM based on the pseudo-likelihood to filter.
In this section, we discuss numerical simulations for an ERGM whose normalization is not known in closed form, which we also apply to real data in Section (ref). We show the concrete applicability of the SD-ERGM approach based on the pseudo-score and its performance as a filter compared with the cross-sectional MCMC estimates of the standard ERGM. The models we consider have two statistics. The first one is the total number of links present in the network. The second statistic is the GWESP, introduced in Section (ref). The ERGM is thus defined by
To test the efficiency of the SD-ERGM, we simulate a known temporal evolution for the parameters and, at each time step, we sample the exact PMF from the resulting ERGMs. Finally, we use the observed change statistics for each time step to estimate two alternative models: a sequence of cross-sectional ERGMs and the SD-ERGM. In what follows, we indicate the values from the DGP of parameter $s$ at time $t$ as $\bar{\theta}_s^{\left(t\right)}$.
We investigate four DGPs similar to those analyzed in Section (ref). We sample and estimate the models 50 times for each DGP. Figure (ref) compares the cross-sectional estimates and the score-driven filtered paths.
Table (ref) reports the RMSE of the GWESP tvps, averaged over the different realizations for the whole sequence $t=1,2,\dots,T$. The SD-ERGM outperforms the cross-sectional ERGM estimates for all the investigated time-varying patterns. Moreover, when the constant DGP is considered, i.e., $\bar{\theta}_1^{\left(t\right)} = \bar{\theta_1} $ and $\bar{\theta}_2^{\left(t\right)} = \bar{\theta_2} $, the average RMSE of the SD-ERGM is larger but comparable, than the correctly specified ERGM that uses all the longitudinal observations to estimate the parameters. The latter result confirms that the SD-ERGM is a reliable and consistent choice even for the static case.
It is worth noticing that, for sampling and cross-sectional inference, we employed the R package ergm that uses state-of-the-art MCMC techniques for both tasks hunter2008ergm. Hence, we compared the SD-ERGM based on the approximate pseudo-likelihood -- both in the definition of the time-varying parameter update and inference of the static parameters -- with a sequence of exact cross-sectional estimates that are in general known to be better performing than the pseudo-likelihood alternative, as mentioned in Section (ref). Even if the cross-sectional estimates are based on the exact likelihood, while the SD approach is based on an approximation, the SD-ERGM remains the best-performing solution. This provides further evidence of the advantages of SD-ERGM as a filtering tool. Finally, the last column of Table (ref) reports the percentage number of times the LM test of calvori2017testing applied to the SD-ERGM correctly classifies the parameters as time-varying (or static for the constant DGP). The test performs correctly in all the cases considered.
To further investigate the proposed SD-ERGM and its version based on the pseudo-likelihood, in this section, we focus on the ERGM having the total number of links and the total number of mutual links as network statistics:
The static version of this model is known as reciprocity $p^\star$ model snijders2002markov. This model is relevant for our discussion because it allows us to compare the SD time-varying extension based on the pseudo-likelihood with the one based on the exact likelihood. Indeed, it is simple enough that the normalizing function is known in closed form, but it has enough structure that its pseudo-likelihood differs from its exact likelihood. The model results in dyads, i.e., pairs of mutual links $(A_{ij}, A_{ji})$, being independent, while the pseudo-likelihood amounts to assuming independent links. Moreover, since its partition function is available in closed form, such a model can be sampled efficiently without resorting to MCMC methods. This allows us to run extensive numerical simulations in reasonable time to investigate the properties of the confidence bands proposed by Ref. \onlinecite{buccheri_etal2018smoother} in the context of SD-ERGM models.
In this section, we will refer to the pseudo-likelihood-based SD-ERGM as PML-SD-ERGM and to the exact likelihood case as ML-SD-ERGM. We compare the capacity of the two models, used as filters, to recover misspecified dynamics using the same approach as in the previous sections, i.e., we simulate a known DGP for $\theta_L^{\left(t\right)}$ and $\theta_M^{\left(t\right)}$. We focus on a DGP where $\theta_L$ and $\theta_M$ follow two independent $AR(1)$ processes, as the one discussed in (ref). Each $AR(1)$ has $\Phi_1 = 0.98$ and $\epsilon \sim N\left(0,\sigma\right)$ with $\sigma = 0.005$. The $\Phi_{0}$ parameters are chosen such that, on average, the network density equals $0.3$, and the fraction of reciprocated pairs is $0.075$. We select this value because it is between the maximum and minimum fraction of reciprocated links possible for a network of density $0.3$, $0$ and $0.3 (N^2-N)/2$, respectively. When comparing results for different network sizes, we keep the density fixed for all network sizes $N$, thus exploring a dense regime\footnote{We found the conclusions of this section to hold also in a sparse network density regime.}. In our numerical experiment, we first sample sequences of synthetically generated observations repeatedly from different specifications of the DGP. We then estimate the PML and ML versions of the SD-ERGM on those observations and filter the tvps. Finally, we quantify their accuracy, with the average RMSE, across 50 samples with respect to the simulated DGP. In Table (ref), we report the RMSE for both PML-SD-ERGM and ML-SD-ERGM, divided by the RMSE of the cross-sectional standard ERGM, for various combinations of network size $N$ and number of observations $T$. It emerges that both versions of SD-ERGM strongly outperform the cross-sectional ERGM. Moreover, the performances of PML-SD-ERGM are similar to the ones of the exact ML-SD-ERGM.
In the final part of this section, we investigate the possibility of using the method of Ref. \onlinecite{buccheri_etal2018smoother}, that we describe in Appendix (ref), to define confidence bands for the parameters filtered with SD-ERGM. The authors characterize the approximation error when the SD approach filters a set of latent parameters whose true DGP is an auto-regressive process. While we refer to the original manuscript for the details, we point out that their procedure rests upon the assumption that the SD filter approximates the true underlying DGP. The authors prove that this approximation becomes exact in the limit of small variance for the latent parameters. Hence, the confidence bands obtained with their method are theoretically guaranteed to be reliable only in this limit. In practice, assessing whether the application of the confidence bands is justified for a given value of the variance of the DGP is appropriate. Numerical experiments can do this to determine their coverage with a simulated DGP. For example, for the model and the DGP considered in this section, we check the coverage of the confidence bands obtained and report the results in Table (ref), for $N=100$.
We find that the coverage of the confidence bands, for both ML-SD-ERGM and PML-SD-ERGM, approaches the nominal value in the limit of large $T$, while for short time series, their coverages are higher than the nominal value. Hence, in small samples, they should be interpreted as having a confidence of at least their nominal values.
After analyzing synthetic data, this section presents two applications to real temporal networks. Our goal is to show the value of SD-ERGM as a methodology to model temporal networks, irrespective of the specific system that a researcher wants to investigate. The two real networks that we consider have been the object of multiple studies in different streams of literature. They have been investigated in the context of ERGMs using different network statistics. We first consider a network of credit relations among Italian banks. The second real-world application focuses on a network of interest for the social and political science community, namely the network of U.S. senators cosponsoring legislative bills.
Our first empirical application is to data from the electronic Market of Interbank Deposit (e-MID). In this market, banks can extend loans to one another for a specified term and/or collateral. Interbank markets are an important point of encounter for banks' supply and demand of extra liquidity. In particular, e-MID has been investigated in many papers iori2008network,finger2013network,mazzarisi2020dynamic,barucca2018organization. Our dataset contains the list of all credit transactions on each day from June 6, 2009, to February 27, 2015. Our analysis investigates the interbank network of overnight loans aggregated weekly. We follow the literature and disregard the size of the exposures, i.e., the weights of the links. We thus consider a link from bank $j$ to bank $i$ present at week $t$ if bank $j$ lent money overnight to bank $i$, at least once during that week, irrespective of the amount lent. This results in a set of $T = 298 $ weekly aggregated networks. For a detailed dataset description, we refer the reader to Ref. \onlinecite{barucca2018organization}.
In recent years, the amount of lending in e-MID has significantly declined. In particular, it abruptly decreased at the beginning of 2012 due to important unconventional measures (Long Term Refinancing Operations) by the European Central Bank that guaranteed an alternative source of liquidity to European banks. The evident non-stationary nature of the evolution of the interbank network is of extreme interest to us. As mentioned in Sections (ref) and (ref), one of the key strengths of SD-ERGM, used as a filter, is precisely the ability to recover such non-stationary dynamics.
In the following, we use the SD beta model for link forecasting. Specifically, we consider the version with a restricted number of static parameters discussed at the end of Sec. (ref). We divide the data set into two samples. We consider rolling windows of $100$ observations and estimate the parameters $\alpha^{\text{out}}$, $\beta^{\text{out}}$, $\alpha^{\text{in}}$ and $\beta^{\text{in}}$ on each one of those rolling windows. For each window, we test the forecasting performances up to $8$ steps ahead (i.e., roughly two months). The forecast works as follows. Assuming that at time $t$, the last date of the rolling window, we have filtered the value for the parameters $ {\overleftarrow{\theta}}^{\left(t\right)} $ and $ {\overrightarrow{\theta}}^{\left(t\right)}$, we plug the estimated static parameters and the matrix ${\mathbf{Y}^{\left(t\right)}}$ in the SD update and compute the tvps $ {\overleftarrow{\theta}}^{\left(t+1\right)} $ and $ {\overrightarrow{\theta}}^{\left(t+1\right)}$. From the latter, we readily obtain the forecast of the adjacency matrix $$ \mathbb{E}\left[{\mathbf{Y}}^{(t+1)}\vert {\overleftarrow{\theta}}^{\left(t+1\right)}, {\overrightarrow{\theta}}^{\left(t+1\right)}\right]\,, $$ where $t+1$ is the first date of the test sample. The $K$-step-ahead forecast for the SD-ERGM model is obtained by simulating the SD dynamics up to $t+K$ $100$ times\footnote{It is worth stressing that the results become stable after $20$ simulations.}, thus obtaining ${\overrightarrow{\theta}}^{(t+K)}_n$ and ${\overleftarrow{\theta}}^{(t+K)}_n$ for $n=1,\dots,100$, and then taking the average of the expected adjacency matrices $\frac{1}{100}\sum_n \mathbb{E}\left[{\mathbf{Y}}^{(t+K)}\vert {\overleftarrow{\theta}}^{(t+K)}_n, {\overrightarrow{\theta}}^{(t+K)}_n\right]$. Given the forecast values, we compute the rate of false positives and false negatives. Then, we drop the first element from the train set and add the first element of the test sample. We repeat the forecasting exercise, estimating the SD-ERGM parameters on the new train set and testing the performance of the new test sample. We name this procedure a rolling estimate and iterate it until the test sample contains the last eight elements of the time series.
Given a forecast for the adjacency matrix, we evaluate the accuracy of the binary classifier by computing the Receiving Operating Characteristic (ROC) curve. All results are collected and presented in Fig. (ref). The left panel reports the ROC curve for one-step-ahead link forecasting obtained according to the SD-ERGM rolling estimate. The panel also shows three other curves based on the static beta model. Specifically, the green curve results from a naive prediction, where a link tomorrow is forecasted, assuming that the $t+1$ ERGM parameter values are equal to those estimated at time $t$. Once the sequence of cross-sectional estimates of the static ERGM is completed, we take the estimated values $ \widehat{{\overleftarrow{\theta}}}^{\left(t\right)}$ and $ \widehat{{\overrightarrow{\theta}}}^{\left(t\right)}$ as observed and model their evolution with an auto-regressive model of order one, AR(1). That amounts to assuming $ \widehat{{\overleftarrow{\theta}}}^{\left(t+1\right)} = c_0 + c_1 \widehat{{\overleftarrow{\theta}}}^{\left(t\right)} + \epsilon^{\left(t\right)}$, where $c_0$ and $c_1$ are the static parameters of the AR(1), and $\epsilon^{\left(t\right)}$ is a sequence of i.i.d. normal random variables with zero mean and variance $\sigma^2$. A similar equation holds for the out-degree parameters. Using the observations from the training sample, we estimate the parameters $c_0$, $c_1$, and $\sigma^2$ and use them for a standard AR(1) forecasting exercise on the test sample. The results correspond to the orange curve. It is important to stress that while the SD-ERGM forecast requires one static and one time-varying estimation on the train set, we must estimate the static parameters for each date in the train sample in the latter procedure.
The left plot of Fig. (ref) shows that the naive one-step-ahead forecast, despite its simplicity, provides a reasonable result. The best performance corresponds, however, to the forecast based on the SD-ERGM. The AR(1) static ERGM improves on the naive forecast and is slightly worse than the SD-ERGM. However, as commented before, it is more computationally intensive. More importantly, the right panel of Fig. (ref) presents a multi-step-head forecasting analysis result. It emerges that the naive forecast's performance (blue curve), tested up to $K=8$, rapidly deteriorates. In contrast, the SD-ERGM multi-step forecast remains the best performing \footnote{In all the results on link forecasting -- one- or multi-step-ahead -- we excluded the links that are always zero, i.e., they never appear in the train and test samples. The reason is that those are extremely easy to predict, and keeping them would give an unrealistically optimistic picture of the predictability of links in the data set. Notably, the ranking of the methods remains unaltered when we keep all links for performance evaluation.}.
Networks describing the U.S. Congress' bills have been the object of multiple studies fowler2006connecting,faust2002comparing,zhang2008community,cranmer2011inferential,moody2013portrait,wilson2019modeling,lee2020varying. It is thus an appropriate real system for our second application of the SD-ERGM framework. In particular, we want to show that the update rule based on the pseudo-score defined in (ref) can be concretely applied to a real network and draws a different picture when compared to the sequence of cross-sectional ERGM estimates. To build the network, we use the freely available data of voting records in the US Senate voteview covering the period from 1867 to 2015, for a total of 74 Congresses. We define the network of co-voting following Refs. \onlinecite{roy2017change} and \onlinecite{lee2020varying}, where a link between two senators indicates that they voted in agreement on over 75% of the votes among those held in a given senate when they were both present. This procedure results in 74 networks, one for each Congress, starting from the $40$th.
We consider the SD-ERGM with the two network statistics discussed in Section (ref) for this empirical application. As defined in (ref), parameter $\theta_1^{\left(t\right)}$ is associated with the number of edges, while $\theta_2^{\left(t\right)}$ with the GWESP statistic. The fact that the number of nodes is not constant over time is not a problem for our application since we do not consider statistics associated with single nodes. That case -- as, for instance, considering the degrees of the beta model -- would require the number of tvps to be different at each time step.
As we did for the numerical simulations and the previous empirical application, we compare our framework with a sequence of standard ERGMs. This empirical exercise does not aim to conclude the specific network at hand. We aim to show that the two approaches return a qualitatively different picture. The choice between the alternative models and combinations of statistics—possibly based on model selection techniques—is beyond the scope of our exercise.
Using the test for temporal heterogeneity based on SD-ERGM, only the parameter $\theta_2$ turns out to be time-varying. Testing the null hypothesis that each parameter is static, we obtain a p-value of $0.1$ for the link density and $10^{-4}$ for GWESP. To check whether the sequence of cross-sectional estimates is consistent with the hypothesis that the parameters remain constant, we estimate the values $\theta_1^c$, $\theta_2^c$ from an ERGM using all observations. This amounts to compute $\theta^c = \underset{\theta}{\operatorname{arg}\,\operatorname{max}}\; \sum_{t=1}^{74} \log P\left({\mathbf{Y}}^{\left(t\right)},\theta\right)$. Then, for each sequence of cross-sectional estimates $\theta_1^{\left(t\right)}$ and $\theta_2^{\left(t\right)}$, we test the hypothesis of them being normally distributed around the constant values with unknown variance. The p-values resulting from the t-tests are $1.4\times 10^{-6}$ and $0.03$ for parameters $\theta_1$ and $\theta_2$, respectively. This simple test confirms that the two approaches imply quantitatively different parameter behaviors. This emerges from Fig. (ref) that reports the estimates from the SD-ERGM (thick red lines), with their respective 95% confidence intervals (shaded red bands), as well as the cross-sectional ERGM estimates -- one per date (blue dots) or using the entire sample (dashed blue line).
To compute the confidence bands as in Ref. \onlinecite{buccheri_etal2018smoother}, we numerically checked whether the data is compatible with a DGP with a small variance. In practice, we first estimate the SD-ERGM. Then we quantify the variance of the latent parameters by estimating an AR(1) on the filtered time series\footnote{We extensively tested via simulation that, for the model at hand and $T$ and $N$ taken from the data, estimating the variance of the latent parameters in such a way results, on average, in a small underestimation. In checking the coverage of the confidence bands, we considered a DGP with variance increased, with respect to the one estimated on the filtered time series, to compensate for this bias.}. Finally, we repeatedly simulate such an $AR(1)$ DGP, similarly to what is done at the end of section (ref), and check the coverage of the confidence bands. We find that, for the current application, the coverage of the confidence bands is 99.9%, hence larger than the nominal value. These simulation-based results support the reliability of the approximate SD filter and provide a conservative estimate of the confidence bands. This allows us, for example, to deduce that the data is not compatible with a model where one of the two parameters is zero.
In this paper, we proposed a framework for describing temporal networks that extends the well-known Exponential Random Graph Models. In the new approach, the parameters of the ERGM have stochastic dynamics driven by the conditional likelihood score. If the latter is unavailable in closed form, we showed how to adapt the score-driven updating rule to a generic ERGM by resorting to the conditional pseudo-likelihood. In this way, our approach can describe the dynamic dependence of the PMF from virtually all the network statistics usually considered in ERGM applications. We investigated two specific ERGM instances using an extensive Monte Carlo analysis of the SD-ERGM reliability as a filter for tvps. The chosen examples allowed us to highlight the applicability of our method to models with a large number of parameters and to models for which the normalization of the PMF is not available in closed form. The numerical simulations proved the clear superior performance of the SD-ERGM over a sequence of standard cross-sectional ERGM estimates. This is true not only in the sparse network regime but also in the dense case when the number of nodes is far from the asymptotic limit. Finally, we run two empirical exercises on real network data. The first application to the e-MID interbank network showed that the SD-ERGM provides a quantifiable advantage in a link forecasting exercise over different time horizons. The second example of the U.S. Congress co-voting political network enlightened that the ERGM and the SD-ERGM could provide a significantly different picture describing the parameter dynamics.
Our work opens several possibilities for future research. First, the applicability of the test for parameter instability in the context of SD-ERGM with multiple network statistics could be investigated much further. This would require an in-depth analysis of the multi-collinearity issues intrinsic to the ERGM context. Second, the SD-ERGM could be applied to multiple instances of real-world temporal networks. An interesting application would be the study of networks describing the dynamical correlation of neural activity in different parts of the brain karahanouglu2017dynamics. In this context, applying the static ERGM has already proven highly successful simpson2011exponential. The last future development we plan to explore is extending the score-driven framework to the description of weighted temporal networks. Regretfully, this setting deserves more attention in the literature giraitis2016estimating. Still, it is of extreme relevance, particularly from the financial stability perspective and its implications for systemic risk.
The authors have no conflicts to disclose.
We are not allowed to share the data from the e-MID inter-bank market, as a confidentiality agreement between the Scuola Normale Superiore of Pisa and the data provider, LIST SpA, binds us. The interested reader can contact \begingroup \catcode `\\12\catcode `\$12\catcode `&12\catcode `\#12\catcode `\^12\catcode `_12\catcode `%12\relax \endgroupLIST SpA for inquiries concerning data access. The second application's data are publicly available at \begingroup\catcode `\\12\catcode `\$12\catcode `&12\catcode `\#12\catcode `\^12\catcode `_12\catcode `%12\relax \endgroup\endgroupURL .