Score-Driven Exponential Random Graphs: A New Class of Time-Varying Parameter Models for Dynamical Networks

TL;DR

This paper introduces score-driven ERGMs, a novel class of models that capture the time-varying nature of network parameters, improving the modeling and prediction of dynamic networks in finance and politics.

Contribution

It extends ERGMs with score-driven dynamics, enabling flexible modeling of evolving network structures and outperforming static models in real-world applications.

Findings

Enhanced prediction of future links in financial networks.

Effective discrimination between static and dynamic network parameters.

Demonstrated flexibility of SD-ERGMs as data-generating processes.

Abstract

Motivated by the increasing abundance of data describing real-world networks that exhibit dynamical features, we propose an extension of the Exponential Random Graph Models (ERGMs) that accommodates the time variation of its parameters. Inspired by the fast-growing literature on Dynamic Conditional Score models, each parameter evolves according to an updating rule driven by the score of the ERGM distribution. We demonstrate the flexibility of score-driven ERGMs (SD-ERGMs) as data-generating processes and filters and show the advantages of the dynamic version over the static one. We discuss two applications to temporal networks from financial and political systems. First, we consider the prediction of future links in the Italian interbank credit network. Second, we show that the SD-ERGM allows discriminating between static or time-varying parameters when used to model the U.S. Congress co-voting network dynamics.

Tables

Table 1. Table 1: RMSEs (on a percentage base) of the filtered paths averaged over all time-varying parameters and all Monte Carlo replicas of the numerical experiment. Left column: results from the cross sectional estimates of the beta model; right column: score-driven beta model results. Each row correspond to one of the four DGPs.

DGPs	Average RMSE
	beta model	SD-beta model
Const	1.75	0.20
Sin	2.76	0.34
Steps	2.46	0.28
AR(1)	1.82	0.24

Table 2. Table 2: First four columns: RMSEs for the filtered paths of the time-varying parameters, averaged over 50 repetitions, for the evolutions of Figure 2 . The last three columns describe the accuracy of the test for dynamics in the parameters, considering the DGPs in Figure 2 , as well as alternative DGPs where only one parameter is time varying. We report the percentage of times that the LM test correctly identifies the parameter as time-varying (or static in the case of the first DGP). The chosen threshold for the p-values is 0.05 0.05 0.05 .

DGP	Average RMSE				LM Test
	ERGM		SD-ERGM		% Correct Results
	${\theta }_1^{\left(t\right)}$	${\theta }_2^{\left(t\right)}$	${\theta }_1^{\left(t\right)}$	${\theta }_2^{\left(t\right)}$	$({\theta }_1^{\left(t\right)},{\theta }_2)$	$({\theta }_1,{\theta }_2^{\left(t\right)})$
Const	0.02	0.1	0.0006	0.004
Sin	0.02	0.04	0.003	0.005	94%	93%
Steps	0.02	0.03	0.01	0.001	92%	96%
AR(1)	0.02	0.2	0.007	0.01	93%	90%

Table 3. Table 3: RMSEs of ML-SD-ERGM and PML-SD-ERGM, relative to that of the cross sectional ERGM. The averages are obtained over 50 repetitions, for the A ​ R ​ ( 1 ) 𝐴 𝑅 1 AR(1) DGP, described in the text. For each value of T 𝑇 T and N 𝑁 N we report the RMSE of the SD-ERGMs divided by the RMSE of the cross sectional ERGM.

	PML-SD-ERGM			ML-SD-ERGM
	50	100	500	50	100	500
100	0.016	0.011	0.006	0.015	0.011	0.006
300	0.015	0.011	0.006	0.014	0.011	0.007
600	0.014	0.012	0.007	0.014	0.011	0.006

Table 4. Table 4: Coverages of the 95% confidence bands averaged over 50 repetitions, for the A ​ R ​ ( 1 ) 𝐴 𝑅 1 AR(1) DGP, described in the text, and N = 100 𝑁 100 N=100 .

$T$	ML-SD-ERGM	PML-SD-ERGM
300	99.1 %	99.9 %
3000	94.5%	95.7 %

Equations

$$\log P\left(\mathbf{Y}\right)=\sum _s{\theta }_s{h}_s\left(\mathbf{Y}\right)-\log \left(\mathcal{K}\left(\theta \right)\right),$$

$$PL\left(\mathbf{Y}\right)=\prod _{ij}{\pi }_{ij}^{Y_{ij}}{\left(1-{\pi }_{ij}\right)}^{\left(1-Y_{ij}\right)}$$

$$f^{\left(t+1\right)}=w+\mathbf{\beta }{f}^{\left(t\right)}+\mathbf{\alpha }{\mathbf{S}}^{\left(\mathbf{t}\right)}{\nabla }^{\left(t\right)},$$

$$\left(\hat{w},\hat{\mathbf{\beta }},\hat{\mathbf{\alpha }}\right)=\underset{\left(w,\mathbf{\beta },\mathbf{\alpha }\right)}{\mathrm{argmax}}\sum _{t=1}^T\log P\left(y^{\left(t\right)}|f^{\left(t\right)}\left(w,\mathbf{\beta },\mathbf{\alpha },{\left\{y^{\left(t^{\prime }\right)}\right\}}_{t^{\prime }=1}^{t-1}\right)\right).$$

$${\nabla }_s^{\left(t\right)}\left(\theta \right)=\frac{\partial \log P\left({\mathbf{Y}}^{\left(t\right)}|\theta \right)}{\partial {\theta }_s}=h_s\left({\mathbf{Y}}^{\left(t\right)}\right)-\frac{\partial \log \mathcal{K}\left(\theta \right)}{\partial {\theta }_s}.$$

$${\theta }^{\left(t+1\right)}=w+\mathbf{\beta }{\theta }^{\left(t\right)}+\mathbf{\alpha }{\mathbf{S}}^{\left(\mathbf{t}\right)}{\nabla }^{\left(t\right)}\left({\theta }^{\left(t\right)}\right).$$

$$\log P\left({\left\{{\mathbf{Y}}^{\left(t\right)}\right\}}_{t=1}^T|w,\beta ,\alpha \right)=\sum _{t=1}^T\log P\left({\mathbf{Y}}^{\left(t\right)}|{\theta }^{\left(t\right)}\left(w,\beta ,\alpha ,{\left\{{\mathbf{Y}}^{\left(t^{\prime }\right)}\right\}}_{t^{\prime }=1}^{t-1}\right)\right).$$

$$\begin{array}[]{ccc}{\overleftarrow{\theta}}_{s}^{\left(t+1\right)}&=&{w}_{s}^{\text{in}}+{\beta}^{\text{in}}_{s}{\overleftarrow{\theta}}_{s}^{\left(t\right)}+{\alpha}^{\text{in}}_{s}\left(\frac{\sum_{i}\left(Y_{is}^{\left(t\right)}-p_{is}^{\left(t\right)}\right)}{\sqrt{\sum_{i}p_{is}^{\left(t\right)}\left(1-p_{is}^{\left(t\right)}\right)}}\right)\\ &&\\ {\overrightarrow{\theta}}_{s}^{\left(t+1\right)}&=&{w}_{s}^{\text{out}}+{\beta}^{\text{out}}_{s}{\overrightarrow{\theta}}_{s}^{\left(t\right)}+{\alpha}^{\text{out}}_{s}\left(\frac{\sum_{i}\left(Y_{si}^{\left(t\right)}-p_{si}^{\left(t\right)}\right)}{\sqrt{\sum_{i}p_{si}^{\left(t\right)}\left(1-p_{si}^{\left(t\right)}\right)}}\right),\\ \end{array}$$

$${\nabla }^{\left(t\right)}\left(\theta \right)=\frac{\partial \log PL\left({\mathbf{Y}}^{\left(t\right)}|\theta \right)}{\partial {\theta }_s^{\left(t\right)\prime }}=\sum _{ij}{\delta }_{ij}^s\left(Y_{ij}^{\left(t\right)}-\frac{1}{1+e^{-{\sum }_l{\theta }_l{\delta }_{ij}^l}}\right),$$

$$\sum _s{\theta }_s{h}_s\left({\mathbf{Y}}^{\left(t\right)}\right)={\theta }_1\sum _{ij}{Y}_{ij}^{\left(t\right)}+{\theta }_2\text{GWESP}\left({\mathbf{Y}}^{\left(t\right)}\right).$$

$$\sum _s{\theta }_s{h}_s\left({\mathbf{Y}}^{\left(t\right)}\right)={\theta }_L\sum _{ij}{Y}_{ij}^{\left(t\right)}+{\theta }_M\sum _{ij}{Y}_{ij}^{\left(t\right)}{Y}_{ji}^{\left(t\right)}.$$

$$\mathbb{E}\left[{\mathbf{Y}}^{(t+1)}|{\overleftarrow{\theta }}^{\left(t+1\right)},{\overrightarrow{\theta }}^{\left(t+1\right)}\right],$$

Full text

Keywords: Exponential Random Graphs; longitudinal networks; temporal networks; econometric network model; time varying parameters; score driven networks.

1 Introduction

A network, or graph111The two names are used interchangeably in this paper., is a useful abstraction for a system composed by a number of single elements that have some pairwise relation among them. The simplified description of social, economic, biological, transportation systems, often very complex in nature, in terms of nodes and links attracted and still attracts an enormous amount of attention, in a number of different streams of literature (Albert and Barabási, 2002; Bullmore and Sporns, 2009; Newman, 2010; Easley et al., 2010; Allen and Babus, 2011). Formally, a graph $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 labelled 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. In general, links may have an orientation and the corresponding adjacency matrix is not symmetric. In this case, network are dubbed directed networks. Moreover, if the elements of the adjacency matrix are allowed to be different from [math] or $1$, one speaks of weighted networks. In the following, we will focus on directed networks, but will not consider the weighted variant.

Often systems that are fruitfully described as networks evolve in time. When pairwise interactions change over time, one usually speaks of temporal networks (Holme and Saramäki, 2012; Craig and Von Peter, 2014). In the time-varying setting, links can last for a finite interval of time, a quantity usually referred to as duration, or be instantaneous. In the latter case, one speaks of interaction or contact links, and different notations can used to describe them (Rossetti and Cazabet, 2018). In this paper, we will focus on the description of temporal networks as sets of links among nodes evolving in discrete time. Then, a dynamical network is a sequence of networks, each one 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}$. In the following, we will present an approach to time-varying networks that is based on two main ingredients: (i) a parametric probabilistic model, according to which one can sample a network realization, i.e. an adjacency matrix; (ii) a simple mechanism to introduce time-variation on the model parameters and, consequently, to induce a dynamics on the network sequence. Concerning the former point, a natural choice is the class of statistical models for networks, known as Exponential Random Graph Models (ERGMs). As far as point (ii) is concerned, a flexible candidate is suggested by the fast growing literature on the Dynamic Conditional Score-driven models (DCSs). The goal of this paper is to present a new class of models for temporal networks and to provide evidence that the novel approach is versatile and effective in capturing time-varying features. To the best of our knowledge, this is the first time the two frameworks are combined to provide a dynamic description of networks.

Our main contribution is an extension of the ERGM framework that allows model parameters to change over time in a score-driven fashion. The result of our efforts is a class of models for time-varying networks where the information encoded in $\mathcal{F}_{t-1}$ is exploited to filter the time-varying parameters $\theta^{\left(t\right)}$ at time $t$. We refer to this class as Score-Driven Exponential Random Graph Models (SD-ERGMs). At this point, it is worth to comment that a generic SD-ERGM can be also used to generate synthetic sequences of graphs, i.e. it can be considered as a data generating process (DGP). However, we are more inclined to interpret it as an effective filter of latent time-varying parameters, regardless of what the true DGP might be.

The rest of the paper is organized as follows. In Section 2 we review the literature on models for temporal networks. In Section 3 we review some key concepts on ERGMs and observation-driven models. In Section 4 we introduce the new class of models and validate it with extensive numerical experiments for three specific instances of the SD-ERGM. Section 5 presents the results from an application to two real temporal networks: The eMID interbank network for liquidity supply and demand and the U.S. Congress co-voting political network. Section 6 draws the relevant conclusions.

2 Literature Review

To the best of our knowledge, an ERGM with score-driven time-varying parameters has never been considered before. Nevertheless, we are by no means the first to discuss models for time-varying networks. For a review of latent space dynamical network models one can refer to (Kim et al., 2018). Moreover, extensions of the ERGM framework for the description of dynamical networks exist in the literature. Two are the main streams. The first one was pioneered by Robins and Pattison (2001) and subsequently discussed in detail in Hanneke et al. (2010) and Cranmer and Desmarais (2011) and is known as TERGM. This approach builds on the ERGM, but allows the network statistics defining 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. Notably Krivitsky and Handcock (2014) generalized TERGM introducing separable specification, known as STERGM, that allows for separate modeling of incidence and duration of ties. A second approach, more related with our work, allows for the parameters of the ERGM to be time-varying. A notable example of this approach is the Varying-Coefficient-ERGM (VCERGM) proposed in Lee et al. (2020). There, the authors combine the varying-coefficient models’ formalism with ERGM to take into account the possibility of the ERGM parameters to be smoothly time-varying. The approach of VCERGM is different from ours in several respects. First, to infer parameter time-variation at time $t$, it uses all the available observations, including those from future times $t^{\prime}>t$. Thus, in time series jargon, it is a smoother and not a filter. Second, as a consequence, it cannot generate sequences of time-evolving networks. At variance, our approach can be used as a filter, as a DGP for time-varying networks, when Monte Carlo scenarios are required, and, following Buccheri et al. (2018), can be extended as a smoother.

In a related, but different, approach Mazzarisi et al. (2020) consider the possibility of a random evolution of node specific parameters. Notably, the latter work accommodates for sender and receiver effect as well as link persistency. The model can also be used to filter the time-varying parameters in a very specific ERGM. Our contribution differs from the latter in flexibility, for the methods used, and for the general scope. In fact, we discuss and test our approach for a generic ERGM. More importantly, the authors of Mazzarisi et al. (2020) consider the random evolution of the parameters to be driven by an exogenous source of randomness. In this respect, following the language of Cox et al. (1981), they consider a class of model known as parameter driven. As stated above, we consider the dynamics of the parameters to be observation driven, i.e. the innovations are generated from observations, as it will be more clear in Section 4.

Finally, it is important to mention that frameworks alternative to latent space models and temporal extension of ERGM for modeling temporal networks have also been considered in the social science literature. Notable examples are the Stochastic Actor Oriented model (SAOM) of Snijders (1996) and the Relational Event Model (REM) of Butts (2008). For an overview of contributions in the social science community, we refer to the literature therein.

3 Introductory Concepts

Given the different origins of the key ideas that we combine, ERGMs are well known in statistics, social sciences, and physics, while DCSs are a recent development in econometrics, we organize this section in two separate subsections. Additionally, for the reader’s convenience, we complement our paper with supplemental information, hereafter SI, that, in its first section includes more details on the two topics.

3.1 ERGM - Exponential Random Graph Models

A statistical model for graphs can be specified providing the probability distribution over the set of possible graphs, i.e. all possible adjacency matrices (see Kolaczyk, 2009, for an introduction and a review of statistical models for networks). If the distribution belongs to the exponential family (Barndorff-Nielsen, 2014), than the model is named ERGM, and its log-likelihood takes the form

[TABLE]

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 literature on ERGMs is extremely vast and still growing (see Schweinberger et al., 2018, for a recent literature review). It is worth to mention that the ERGM framework is intrinsically linked to the very well known principle of maximum entropy (Shannon, 2001) and its applications to statistical physics (Jaynes, 1957). 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)$ (Park and Newman, 2004; Garlaschelli and Loffredo, 2008). The sufficient statistics, known as network statistics, are functions of the adjacency matrix $h_{s}\left({\mathbf{Y}}\right)$ and the probability mass function (PMF) is defined by (1). The normalizing factor $\mathcal{K}\left(\theta\right)$ is often not available as a closed-form function of the parameters $\theta$. Each matrix element is a binary random variable, and its probability depends only on the value of the network statistics appearing in (1).

In the following we will focus on two specific examples of ERGMs that describe distinct features of the network and require different approaches to the parameter inference. The first one is meant to capture the heterogeneity in the number of connections that each node can have and it allows for straightforward maximum likelihood estimation, and is known as beta model, fitness model, and configuration model (Zermelo, 1929; Holland and Leinhardt, 1981; Caldarelli et al., 2002; Park and Newman, 2004; Garlaschelli and Loffredo, 2008; Chatterjee et al., 2011). It is defined as an ERGM having as network statistics the number of connections each node has - its degree. For the directed network case considered here, we have – for node $i$ – out-degree $\overrightarrow{D}_{i}=\sum_{j}{Y_{ij}}$ and in-degree $\overleftarrow{D}_{i}=\sum_{j}Y_{ji}$. Hence this model has two parameters per each node: ${\overrightarrow{\theta}}_{i}$, that captures the propensity of node $i$ to form outgoing connections, and ${\overleftarrow{\theta}}_{i}$ those incoming.

The second one, that we use in Section 4.2, is an ERGM having as statistic the Geometrically Weighted Edgewise Shared Partners (GWESP). That is a network statistic that describes 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 Snijders et al. (2006); Robins et al. (2007) and discussed in Hunter and Handcock (2006).

While for the beta model standard maximum likelihood estimation (MLE) can be carried out (Chatterjee et al., 2011), often inference in ERGMs is complicated by the fact that the normalizing factor in (1) as a function of the parameters is not available in closed-form. As is the case for an ERGM including GWESP as network statistic. In such cases, there are two standard approaches to ERGM inference, both consisting in maximizing alternative functions that are known to share the same optimum as the exact likelihood. The first possibility (described, for example, in Snijders, 2002) is to maximize an objective function obtained from a sufficiently large sample drawn from the PMF with an arbitrary (but close enough to the true one) parameter. As a consequence of the non independence of the links in the general ERGM, sampling from (1) necessary relies on Markov Chain Monte Carlo (MCMC) approaches (see Hunter et al., 2008, for a description of a popular software that implements it). The computational burden of MCMC-based estimation can be prohibitive for large enough graphs. For this reason, a second approximate inference procedure, known as Maximum Pseudo-Likelihood Estimation (MPLE), first proposed for ERGMs in the seminal work of Strauss and Ikeda (1990), is often used in empirical applications. MPLE is based on the optimization of the pseudo-likelihood function, that is in turn defined from link specific variables (one for each element of the adjacency matrix) named change statistics. Given an ERGM, the change statistic for the link between node $j$ and $i$, associated with network statistic $h_{s}$ is $\delta^{s}_{ij}=h_{s}\left({\mathbf{Y}}_{ij}^{+}\right)-h_{s}\left({\mathbf{Y}}_{ij}^{-}\right)$, where ${\mathbf{Y}}_{ij}^{+}$ is a matrix such that $Y^{+}_{ij}=1$ and it is equal to ${\mathbf{Y}}$ in all other elements. Similarly, ${\mathbf{Y}}_{ij}^{-}$ has $Y^{-}_{ij}=0$ and it is equal to ${\mathbf{Y}}$ in all other entries. Given these definitions, the pseudo-likelihood reads

[TABLE]

where $\pi_{ij}=\left(1+e^{-\sum_{s}\theta_{s}\delta^{s}_{ij}}\right)^{-1}$. In Section $1.2$ of the SI we give more details on GWESP and ERGM inference based on the pseudo-likelihood.

3.2 Score-Driven Models

The second main ingredient of this work is the class of Dynamic Conditional Score-driven models introduced in Creal et al. (2013) and Harvey (2013), also known as Generalized Autoregressive Score (GAS) models 222see http://www.gasmodel.com/index.htm for the updated collection of papers dealing with GAS models.. In the language of Cox et al. (1981), DCSs belong to the class of observation-driven models. Specifically, in order to introduce DCSs, 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)}|f^{\left(t\right)}\right)$, that depends on a vector of time-varying parameters $f^{\left(t\right)}\in\mathbb{R}^{K}$. Defining the score as $\nabla^{\left(t\right)}=\frac{\partial\log{P\left(y^{\left(t\right)}|f^{\left(t\right)}\right)}}{\partial f^{\left(t\right)\prime}}$, a score-driven model assumes that the time evolution of $f^{\left(t\right)}$ is ruled by the recursive relation

[TABLE]

where $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, that is often chosen to be the inverse of the square root of the Fisher information matrix associated with $P\left(y^{\left(t\right)}|f^{\left(t\right)}\right)$.However, this is not the only possible specification and different choices for the scaling are discussed in Creal et al. (2013).

The most important feature of (3) is the role of the score as a driver of the dynamics of $f^{\left(t\right)}$. 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 a stochastic dynamics exactly thanks to the random occurrence of $y^{\left(t\right)}$. When the score-driven recursion is regarded as a filter, the update rule in (3) is used to obtain a sequence of filtered $\left\{\hat{f}^{\left(t\right)}\right\}_{t=1}^{T}$. In this setting, the static parameters are estimated maximizing the log-likelihood of the whole sequence of observations (for a detailed discussion, see Harvey, 2013; Blasques et al., 2014),

[TABLE]

A second look at eq. (3) 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. Indeed, at each step the score pushes the parameter vector along the log-likelihood steepest direction. After scaling with the matrix $\boldsymbol{S}$, the intensity of the push is modulated by the parameter $\boldsymbol{\alpha}$, and its direction adjusted by the auto-regressive component. Moreover, there are motivations, originating in information theory, for the optimality of the score-driven updating rule Blasques et al. (2015), as we discuss in Section $2$ of the SI .

It turns out that many well known models in econometrics can be expressed as score-driven models. Famous examples are the Generalized Autoregressive Conditional Heteroskedasticity (GARCH) model of Bollerslev (1986), the Exponential GARCH model of Nelson (1991), the Autoregressive Conditional Duration (ACD) model of Engle and Russell (1998), and the Multiplicative Error Model (MEM) of Engle (2002). The introduction of this framework in its full generality opened the way to applications in various contexts.

In order to provide confidence bands for the filtered parameters, Blasques et al. (2016) discussed methods to quantify the uncertainty associated with the score driven filters, when the DGP is itself a score driven model. Their method is meant to capture only the uncertainty due to the estimation of the static parameters, often referred to as parameter uncertainty. Hence the confidence bands reliably quantify uncertainty only when the DGP is score driven. In other words, these bands do not take into account what is known as filtering uncertainty. This is the uncertainty due to the fact that in general we do not know the true DGP and the score-driven filter may be regarded only as an approximate filter. Recently Buccheri et al. (2018) investigated the approximation error made by applying a score driven filter to a time varying parameter model following a different DGP. They found that, for a class of DGPs where the parameters follow an auto-regressive process, the approximation becomes exact in the limit of small variance of the latent parameters. Moreover they proposed a method to define confidence bands, inspired by Hamilton (1986), that accounts for both filtering and parameter uncertainty in Score Driven filters.

As a final aspect, which will be relevant in the application section of this paper, score-driven models allow for a test discriminating whether the observations are better described by a model with time-varying parameters or static ones. In fact, following Engle (1982), Calvori et al. (2017) discuss the performances of a test for parameter temporal variation tailored for score-driven models, of which we give some more details in Section $2$ of the SI.

4 Score-Driven Exponential Random Graphs

In this section, we introduce the general SD-ERGM framework, discuss in detail the applicability of the score-driven approach to three different ERGMs, and validate their performances with extensive numerical simulations.

We propose to apply the score-driven methodology to ERGMs, in order to allow any of the parameters $\theta_{s}$ in (1) 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 the description of time-varying networks, more than in a single model, in very much the same way as ERGM is considered a modeling framework for static networks. We refer to such class of models as Score-Driven Exponential Random Graphs Models.

Conceptually, applying the score-driven approach is fairly straightforward. Given the observations $\left\{{Y_{ij}^{\left(t\right)}}\right\}_{t=1}^{T}$, we can apply the update rule in (3) to all or some elements of $\theta$, each of which is associated with a network statistic in (1). In order 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

[TABLE]

It follows that the vector of time-varying parameters evolves according to

[TABLE]

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 a stochastic dynamics because, at each time $t$, the adjacency matrix is not known in advance. It is randomly sampled from $P\left({\mathbf{Y}^{\left(t\right)}}|\theta^{\left(t\right)}\right)$ and then used to compute the score that, as a consequence, becomes itself stochastic. 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 estimated maximizing the log-likelihood of the whole time series. Taking into account that each network ${\mathbf{Y}^{\left(t\right)}}$ is independent from all the others conditionally on the value of $\theta^{\left(t\right)}$, the log-likelihood can be written as

[TABLE]

It is evident that the computation of the normalizing factor, and its derivative with respect to the parameters, is essential for the SD-ERGM. Not only it enters the definition of the update, but it is also required for the optimization of (6).

Our main motivation for the introduction of SD-ERGM is to describe the time evolution of a sequence of networks by means of the evolution of the parameters of an ERGM. We assume to know, from the context or from previous studies of static networks in terms of ERGM, which statistics are more appropriate in the description of a given network. Hence, we do not discuss the choice of statistics in the context of dynamical networks, but refer the reader to Goodreau (2007) and Hunter et al. (2008) for examples of feature selection and Goodness Of Fit (GOF) evaluation, as well as to Shore and Lubin (2015) for a recent proposal to quantify GOF specifically in network models.

In the rest of this section, we discuss in details the SD extension of ERGMs with given statistics. The first example allows for the exact computation of the likelihood, but the number of parameters can become large for large network. In the second example, we discuss 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 time-varying parameters when the DGP is known and the score-driven model is employed as a misspecified filter. Moreover, we show a first application of the LM test in assessing the time-variation of ERGM parameters.

4.1 Score Driven Beta Model

Our first specific example is the score-driven version of the beta model, introduced in Sec. 3.1 and discussed in detail in Section $1.1$ of the SI. We start with this model not only because of its wide applications and relevance in various streams of literature, but also 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 need to describe the dynamics of a large amount of parameters, this last feature poses a challenge to 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 a high dimensional dynamics.

As anticipated, the SD-beta model is defined by the application of (5) to each one 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 very common, in score-driven models with numerous time-varying parameters, to restrict the matrices $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ of (3) 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

[TABLE]

where the superscripts $in$ and $out$ indicate the first and second half of the parameter vectors, respectively. In order to simplify the inference procedure, we consider a two-step approach. First, we fix the node specific parameters $w_{i}$ in order 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.

4.1.1 SD-ERGMs as filters: Numerical Simulations

As mentioned in the Introduction, SD-ERGMs (as other observation driven model, e.g. GARCH) can be seen either as DGP and estimated on real time series or as predictive filters (Nelson (1996)), since time-varying parameters follow one-step-ahead predictable processes. In this Section we show the power of the ERGMs in this second setting. Specifically, we simulate generic non-stationary temporal evolution for the parameters $\theta^{\left(t\right)}$ of temporal networks. We then use the SD-ERGM to filter the paths of the parameters and evaluate its performances. It is important to note that the simulated dynamics of the parameters is different from the score-driven one used in the estimation.

In practice, at each time $t$, we sample the adjacency matrix from the PMF of an ERGM with parameters333In 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 in order to ensure heterogeneity in the expected degrees of the nodes under the static beta model. For the remaining three DGPs, half of the parameters is static and half is time-varying evolving with either a deterministic sinusoidal function, a deterministic step function and a stochastic AR(1) dynamics. More details on the definition of such DGPs are given in Section $3$ of the SI.

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 time-varying parameters and $100$ simulations, and report the results in Table 1. These results confirm that the SD beta model outperforms the standard beta model in recovering the true time-varying pattern. Notably, this holds true even when the DGP is inherently non stationary, as in the case of the DGP where each parameter has a step like evolution. Indeed the results of this Section and of Section 4.2 confirm that, while the SD update rule (3) defines a stationary DGP (see Creal et al., 2013), using SD models as filters, we can effectively recover non stationary parameters’ dynamics.

Our last numerical simulations for the SD beta model explore its applicability and performances for networks of increasing size. We explore this setting for two reasons. The first one is that many real systems are described by networks with a large number of nodes. The second reason is that we want to compare the performances of our approach with those of the standard beta model in regimes where the latter is known to perform better, under suitable conditions. Indeed, as mentioned in Section $1.1$ of the SI, asymptotic results on the single observation estimates Chatterjee et al. (2011) guarantee that, if the network density remains constant as $N$ grows larger, 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. In order 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 constant the average degree. Moreover, to ease the computational burden for the estimates, we consider a restricted version of the SD-Beta model, as detailed in Section $3.1$ of the SI, 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.

In this analysis, we consider 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 we let it depend on $N$ in two distinct ways, each one corresponding to a different density regime: one generating sparse networks and the other dense ones. It is important to notice that the asymptotic results of (Chatterjee et al., 2011) 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 1. Then, for both regimes and each value of $N$, we compute the average RMSE across all time-varying parameters and all Monte Carlo replicas. In the right panel of Figure 1, the average RMSEs for different values of $N$ clearly indicate that, also for large networks, the SD version of the beta model attains better results compared with the cross sectional estimates. As expected, in the dense network regimes, both approaches reach the same accuracy 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.

4.2 Pseudo Likelihood SD-ERGM

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 (3) 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. 3.1, that we refer to as pseudo-score

[TABLE]

in place of the exact score in the definition the SD-ERGM update (3). Additionally, we use the pseudo-likelihood for each observation ${\mathbf{Y}}^{\left(t\right)}$ in (6) for the inference of 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)$ 444For 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 Hunter et al. (2008).. In the following, we show that the update based on the score of the pseudo-likelihood is effective in filtering the path of time-varying parameters. Remarkably, this is true even when the probability distribution in the DGP is the exact one, i.e. when we sample from the exact likelihood and then use the SD-ERGM based on the pseudo-likelihood to filter.

4.2.1 SD-ERGM for Transitivity and Network Density

In this section we discuss numerical simulations for an ERGM whose normalization is not known in closed form, that we apply also to real data in Section 5.2. We show the concrete applicability of the SD-ERGM approach based on the pseudo-score and its performance as a filter in comparison 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 statistics is the GWESP, introduced in Section 3.1. The ERGM is thus defined by

[TABLE]

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 4.1. We sample and estimate the models 50 times for each DGP. Figure 2 compares the cross sectional estimates and the score-driven filtered paths.

Table 2 reports the RMSE of the GWESP time-varying parameters, averaged over the different realizations for the whole sequence $t=1,2,\dots,T$. It is evident that 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, even for the static case, the SD-ERGM is a reliable and consistent choice.

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 (see Hunter et al., 2008, for a description of the software). 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 3.1. 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. In our opinion, this provides further evidence of the advantages of SD-ERGM as a filtering tool. Finally, the last column of Table 2 reports the percentage number of times the LM test of Calvori et al. (2017) 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.

4.2.2 Comparison of Pseudo and Exact Likelihood SD-ERGM

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:

[TABLE]

The static version of this model is known as reciprocity $p^{\star}$ model (Snijders, 2002). 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 such that its pseudo-likelihood differs from its exact likelihood. In fact, 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 the need to resort to MCMC methods. This allows us to run extensive numerical simulations, in reasonable time, to investigate the properties of the confidence bands proposed by Buccheri et al. (2018) 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 a 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 4.1. 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 is equal to $0.3$ and the fraction of reciprocated pairs is $0.075$. We choose this value because it is between the maximum and minimum fraction of reciprocated links possible for a network of density $0.3$, which is [math] and $0.3(N^{2}-N)/2$. When comparing results for different network sizes, we keep the density fixed for all network sizes $N$, thus exploring a dense regime555We found the conclusions of this section to hold true also in a sparse network density regime.. In our numerical experiment, we first sample repeatedly sequences of synthetically generated observations from different specifications of the DGP. We then estimate the PML and ML versions of the SD-ERGM on those observations, and filter the time varying parameters. Finally we quantify their accuracy, with the average RMSE, across 50 samples, with respect to the simulated DGP. In Table 3 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 clearly emerges that both versions of SD-ERGM strongly outperform the cross sectional ERGM. Moreover, the performances of PML-SD-ERGM is similar to the one of the exact ML-SD-ERGM.

In the final part of this section, we investigate the possibility of using the method of Buccheri et al. (2018), that we describe in Section $2$ of the SI, to define confidence bands for the parameters filtered with SD-ERGM. In Buccheri et al. (2018) the authors characterize the approximation error when the SD approach is used to filter a set of latent parameters, whose true DGP is an auto-regressive process. While we refer to Buccheri et al. (2018) 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, it is appropriate to assess whether for a given value of the variance of the DGP the application of the confidence bands is justified. This can be done by numerical experiments to assess 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 4, for $N=100$.

We find that 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 is higher than the nominal value. Hence in small samples they should be interpreted as having a confidence of at least their nominal values.

5 Applications to Real Data

After the analysis of synthetic data, this section presents two applications to real dynamical 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.

5.1 Link Prediction in Interbank Networks

Our first empirical application is to data from the electronic Market of Interbank Deposit (e-MID), a market where 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 (see, for example Iori et al., 2008; Finger et al., 2013; Mazzarisi et al., 2020; Barucca and Lillo, 2018, and references therein). Our dataset contains the list of all credit transactions in each day from June 6, 2009 to February 27, 2015. In our analysis, we investigate 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 description of the dataset, we refer the reader to Barucca and Lillo (2018).

In recent years, the amount of lending in e-MID has significantly declined. In particular, as discussed in Barucca and Lillo (2018), it abruptly declined at the beginning of 2012, as a consequence of 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 inter-bank network is of extreme interest for our purposes. In fact, as mentioned in Sections 4.1 and 4.2, one of the key strengths of SD-ERGM, used as a filter, is precisely the ability of recovering such non-stationary dynamics.

In the following, we use the SD beta model for links forecasting. Specifically, we consider the version with a restricted number of static parameters discussed at the end of Sec. 4.1.1. We divide the data set in 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 then 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 time-varying parameters ${\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

[TABLE]

where $t+1$ is the first date of the test sample. The $K$-step-ahead forecast for the SD-ERGM model is obtained simulating the SD dynamics up to $t+K$ $100$ times666It is worth to stress 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)}|{\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 to it 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 on the new test sample. We name this procedure rolling estimate and iterate it until the test sample contains the last 8 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 Figure 3. 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 other three curves based on the static beta model. Specifically, the green curve results from a naive prediction, where the presence of 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, in the latter procedure we have to estimate the static parameters for each date in the train sample.

The left plot of Fig. 3 shows that the naive one-step-ahead forecast, in spite of its simplicity, provides a quite 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 it is slightly worst than the SD-ERGM. However, as commented before, it is more computationally intensive. More importantly, the right panel of Fig. 3 presents the result from a multi-step-head forecasting analysis. It emerges clearly that the performance of the naive forecast (blue curve), tested up to $K=8$, rapidly deteriorates, while the SD-ERGM multi-step forecast remains the best performing. 777In 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 on the predictability of links in the data set. Importantly, the ranking of the methods remain unaltered when we keep all links for performance evaluation..

5.2 Temporal Heterogeneity in U.S. Congress Co-Voting Political Network

Networks describing the U.S. congress’ bills have been the object of multiple studies (see, for example Fowler, 2006; Faust and Skvoretz, 2002; Zhang et al., 2008; Cranmer and Desmarais, 2011; Moody and Mucha, 2013; Wilson et al., 2016; Lee et al., 2020). 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 (8) can be concretely applied to a real network, and that it draws a different picture when compared to the sequence of cross sectional ERGM estimates. In order to build the network, we use the freely available data of voting records in the US Senate (see Lewis et al., 2019) covering the period from 1867 to 2015, for a total of 74 Congresses. We define the network of co-voting following Roy et al. (2017) and Lee et al. (2020), 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 a sequence of 74 networks, one for each different Congress starting from the $40$th.

For this empirical application, we consider the SD-ERGM with the two network statistics discussed in Section 4.2. As defined in (9), 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 to single nodes. That case – as for instance considering the degrees of the beta model – would require the number of time-varying parameters 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. The goal of this empirical exercise is not to draw conclusions about the specific network at hand. Our main aim is 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. In fact, 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. In order 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 behaviors for the parameters. This clearly emerges from Figure 4 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).

In order to compute the confidence bands as in Buccheri et al. (2018), we checked numerically whether the data is compatible with a DGP having 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 series888We 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 done at the end of section 4.2.2, 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 greater 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. Thus allowing us, for example, to deduce that the data is not compatible with a model where one of the two parameters is zero.

6 Conclusions

In this paper, we proposed a framework for the description of temporal networks that extends the well known Exponential Random Graph Models. In the new approach, the parameters of the ERGM have a stochastic dynamics driven by the score of the conditional likelihood. If the latter is not available 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 by means of an extensive Monte Carlo analysis of the SD-ERGM reliability as a filter for time-varying parameters. 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 not only true 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 networks data. The first application to 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 on the U.S. Congress co-voting political network enlightened that the ERGM and the SD-ERGM could provide a significantly different picture in describing the parameter dynamics.

Our work opens a number of 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 that are intrinsic to the ERGM context. Second, the SD-ERGM could be applied on multiple instances of real world dynamical networks. An interesting application would be the study of networks describing the dynamical correlation of neural activity in different parts of the brain (see, for example, Karahanoğlu and Van De Ville, 2017, for a review of the topic and list of references). In this context, the application of the static ERGM have already proven to be extremely successful (as, for example, in Simpson et al., 2011). The last future development that we plan to explore is the extension of the score-driven framework to the description of weighted dynamical network. Regretfully, this setting has not received enough attention in the literature (one isolate example is Giraitis et al., 2016), but it is of extreme relevance, particularly from the financial stability perspective and its implications for systemic risk.