Regression-based Network Reconstruction with Nodal and Dyadic Covariates and Random Effects

by   Michael Lebacher, et al.

Network (or matrix) reconstruction is a general problem which occurs if the margins of a matrix are given and the matrix entries need to be predicted. In this paper we show that the predictions obtained from the iterative proportional fitting procedure (IPFP) or equivalently maximum entropy (ME) can be obtained by restricted maximum likelihood estimation relying on augmented Lagrangian optimization. Based on the equivalence we extend the framework of network reconstruction towards regression by allowing for exogenous covariates and random heterogeneity effects. The proposed estimation approach is compared with different competing methods for network reconstruction and matrix estimation. Exemplary, we apply the approach to interbank lending data, provided by the Bank for International Settlement (BIS). This dataset provides full knowledge of the real network and is therefore suitable to evaluate the predictions of our approach. It is shown that the inclusion of exogenous information allows for superior predictions in terms of L_1 and L_2 errors. Additionally, the approach allows to obtain prediction intervals via bootstrap that can be used to quantify the uncertainty attached to the predictions.



There are no comments yet.


page 1

page 2

page 3

page 4


Iterative Estimation of Mixed Exponential Random Graph Models with Nodal Random Effects

The presence of unobserved node specific heterogeneity in Exponential Ra...

Estimating the size of a closed population by modeling latent and observed heterogeneity

The paper describe a new class of capture-recapture models for closed po...

Median bias reduction in random-effects meta-analysis and meta-regression

Random-effects models are frequently used to synthesise information from...

On the benefits of maximum likelihood estimation for Regression and Forecasting

We advocate for a practical Maximum Likelihood Estimation (MLE) approach...

Separable correlation and maximum likelihood

We consider estimation of the covariance matrix of a multivariate normal...

Estimation of Latent Network Flows in Bike-Sharing Systems

Estimation of latent network flows is a common problem in statistical ne...

Incorporating Actor Heterogeneity into Large Network Models through Variational Approximations

The analysis of network data has gained considerable interest in the rec...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The problem of how to obtain predictions for unknown entries of a matrix, given restrictions on the row and column sums is a problem that comes with many labels. Without a sharp distinction of names and fields, some non exhaustive examples for keywords that are related to very similar settings are Network Tomography and Traffic Matrix Estimation

in Computer Sciences and Machine Learning (e.g. 

Vardi, 1996, Coates et al., 2002, Hazelton, 2010, Airoldi and Blocker, 2013, Zhou et al., 2016), Input-Output Analysis in Economics (e.g. Bacharach, 1965, Miller and Blair, 2009), Network Reconstruction in Finance and Physics (e.g. Sheldon and Maurer, 1998, Squartini and Garlaschelli, 2011, Mastrandrea et al., 2014, Gandy and Veraart, 2018), Ecological Inference in Political Sciences (e.g. King, 2013, Klima et al., 2016), Matrix Balancing in Operation Research (e.g. Schneider and Zenios, 1990) and many more.

An old but nevertheless popular solution to problems of this kind is the so called iterative proportional fitting procedure (IPFP), firstly introduced by Deming and Stephan (1940) as a mean to obtain consistency between sampled data and population-level information. In essence, this simple procedure iteratively adjusts the estimated entries until the row and column sums of the estimates match the desired ones. In the statistics literature, this procedure is frequently used as a tool to obtain maximum likelihood estimates for log-linear models in problems involving three-way and higher order tables (Fienberg et al., 1970, Bishop et al., 1975, Haberman, 1978, 1979). Somewhat parallel, the empirical economics literature, concerned with the estimation of Input-Output matrices, proposed a very similar approach (Bacharach, 1970), often called RAS algorithm. Here, the entries of the matrix must be consistent with the inputs and the outputs. The solution to the problem builds on the existence of a prior matrix that is iteratively transformed to a final matrix that is similar to the initial one but matches the input-output requirements. Although the intention is somewhat different, the algorithm is effectively identical to IPFP (Onuki, 2013). The popularity of the procedure can also be explained by the fact that it provides a solution for the so called maximum entropy (ME) problem (Malvestuto, 1989, Upper, 2011, Elsinger et al., 2013). In Computer Sciences, flows within router networks are often estimated using Raking and so called Gravity Models (see Zhang et al., 2003). Raking is in fact identical to IPFP and the latter can be interpreted as a special case of the former.

In this paper, we propose an estimation algorithm that builds on augmented Lagrangian optimization (Powell, 1969, Hestenes, 1969) and can provide the same predictions as IPFP but is flexible enough to be extended toward more general concepts. In particular we propose to include exogenous covariates and random effects to improve the predictions of the missing matrix entries. Furthermore, we compare our approach with competing models using real data. To do so, we look at an international financial network of claims and liabilities where we pretend that the inner part of the matrix is unknown. Since in the data at hand the full matrix is in fact available we can carry out a competitive comparison with alternative routines. Note that commonly the inner part of the financial network remains unknown but finding good estimates for the matrix entries is essential for central banks and financial regulators. This is because it is a necessary prerequisite for evaluating systemic risk within the international banking system. See e.g. a very recent study by researchers from 15 different central banks (Anand et al., 2018) where the question of how to estimate financial network linkages was identified as being crucial for contagion models and stress tests. Our proposal has therefore a direct practical contribution.

The paper is structured as follows. In Section 2 we relate maximum entropy, maximum likelihood and IPFP. In Section 3 we introduce our model and discuss estimation and inference as well as potential extensions. After a short data description in Section 4 we apply the approach, compare different models and give a small simulation study that shows properties of the estimator. Section 5 concludes our paper.

2 Modelling approach

2.1 Notation

Our interest is in predicting non-negative, directed dyadic variables among observational units at time points . The restriction to non-negative entries is henceforth referred to as non-negativity constraint. We do not allow for self-loops and leave elements undefined. Hence, the number of unknown variables at each time point is given by . Let be an

-dimensional column vector and define

as the corresponding ordered index set. We denote the th row sum by and the th column sum by . Stacking the row and column sums, results in the -dimensional column vector . Furthermore, we define the known binary routing matrix such that the linear relation


holds. Henceforth, we will refer to relation (1) as marginal restrictions. Furthermore, we denote each row of by the row vector . Hence, we can represent the marginal restrictions row wise by

Note that in cases where some elements of are zero, the number of unknown variables to predict decreases and matrix

must be rearranged accordingly. In the following we will ignore this issue and suppress the time-superscript for ease of notation. Random variables and vectors are indicated by upper case letters, realizations as lower case letters.

2.2 Maximum entropy, iterative proportional fitting and maximum likelihood

Besides the long known relation between maximum entropy (ME) and IPFP, there also exists an intimate relation between maximum entropy and maximum likelihood that is formalized for example by Golan and Judge (1996) and is known as the Duality theorem, see for example Brown (1986) and Dudík et al. (2007). Also in so-called configuration models (Squartini and Garlaschelli, 2011, Mastrandrea et al., 2014) the connection between maximum entropy and maximum likelihood is a central ingredient for network reconstruction.

In the following (i) we rely on the work of Golan and Judge (1996), Squartini and Garlaschelli (2011) and Muñoz-Cobo et al. (2017)

in order to briefly derive the ME-distribution in the given setting. (ii) After that, we show that IPFP indeed maximizes the ME-distribution. (iii) Based on the first two results, we show that we can arrive at the same result as IPFP by constrained maximization of a likelihood where each matrix entry comes from an exponential distribution.

(i) Maximum entropy distribution: We formalize the problem by defining the Shannon entropy functional of the system as

where we make it explicit in the notation that the functional takes the function as input. The support of is given by , ensuring the non-negativity constraint. Furthermore, we require that the density function integrates to unity


We denote the expectation of the random vector by and formulate the marginal restrictions in terms of linear restrictions on which we specify as


Combining the constraints (2) and (3) results into the Lagrangian functional


with Lagrange multipliers for . The solution can be found using the Euler-Lagrange equation (Dym and Shames, 2013), stating that a functional of the form is stationary (i.e. its first order derivative is zero) if


If the Lagrangian functional does not depend on the derivative of , we find the right hand side in equation (5) to be zero so that no derivative appears. For the Lagrangian functional (4) this provides


Rearranging the terms in (6) results in the maximum entropy distribution


In order to ensure restriction (2) we set where is the parameter vector and

where for ensures integration to a finite value. Taken together, this leads to the exponential family distribution


Apparently, the sufficient statistics in (8) result through

and hence, we can characterize the dimensional random variable in terms of parameters . Using (8), the second order condition results from

and ensures that is indeed a maximizer.

(ii) IPFP and the maximum entropy distribution: In order to solve for the parameters of the maximum entropy distribution we take the first derivative of the log-likelihood obtained from (8), i.e.


Since (8) is an exponential family distribution we can use the relation

and the maximum likelihood estimator results from the score equations


If we now multiply the left and the right hand side of (10) by parameter we get

and can solve the problem using fixed-point iteration (Dahmen and Reusken, 2006). That is we fix the right hand side to and update the left side to through


But this is in fact iterative proportional fitting, a procedure that iteratively rescales the parameters until the estimates match the marginal constraints. Convergence is achieved when , satisfying the score equations (10). More generally, the log-likelihood (9) is monotonically non-decreasing in each update step (11) and convergence of (11) is achieved only if the log-likelihood is maximized (Koller et al., 2009, Theorem 20.5).

(iii) IPFP and constrained maximum likelihood: If we re-sort the sufficient statistics and re-label the elements of we get

with for . This leads to


and where with (10) . Hence, we can represent the whole system as the product of densities from exponentially distributed random variables for . That is


with observed margins and .

3 Maximum likelihood-based estimation strategy

3.1 Parametrization and estimation

From result (13) it follows that we can use a distributional framework in order to build a generalized regression model. We exemplify this with a model which includes a sender-effect denoted as and a receiver-effect . We stack the coefficients in a parameter vector such that the following log-linear expectation results


Based on this structural assumption, we can now maximize the likelihood derived from (13) with respect to subject to the observed values and the moment condition

. In the given formulation, the moment condition is linear in

but not in . Consequently, the numerical solution to the problem might be burdensome. We therefore propose to use an iterative procedure that is somewhat similar to the Expectation Conditional Maximization (ECM, Meng and Rubin, 1993) algorithm, since it involves iteratively forming the expectation of based on the previous parameter estimate (E-Step) and constrained maximization afterwards (M-Step). To be specific, define the design matrix , that contains indices for sender- and receiver-effects. Matrix has rows indexed by where for we have the -th and the-th element of equal to and all other elements are equal to zero. Starting with an initial estimate that satisfies , we form the expectation of the log-likelihood

Then, the maximization problem in the M-step is given by


A suitable optimizer for non-linear constraints is available by the augmented Lagrangian (Hestenes, 1969, Powell, 1969)


with and being auxiliary parameters. The augmented Lagrangian method decomposes the constrained problem (15) into iteratively solving unconstrained problems. In each iteration we start with an initial parameter in order to find the preliminary solution . Then, we update in order to increase the accuracy of the estimate. In case of slow convergence, also can be increased. An implementation in R is given by the package nloptr by Johnson (2014).

3.2 Confidence and prediction intervals

Considering the data entries as exponentially distributed allows for a quantification of the uncertainty of the estimates. We pursue this by bootstrapping (Efron and Tibshirani, 1994) here. Given a converged estimator we can draw for each matrix entry from an exponential distribution with expectation in order to obtain bootstrap samples . For each bootstrap sample we calculate the marginals and re-run the constrained estimation procedure resulting in vectors of estimated means . Consequently, the moment condition holds for all mean estimates of the bootstrap and by model-construction, the expected marginal restrictions from the bootstrap sample match the observed ones:

Based on the bootstrap estimates, we can easily derive confidence intervals for

using the variability of for . Additionally, we define the prediction error as and construct prediction intervals for the unknown

based on the quantiles of the empirical distribution of

3.3 Extensions with exogenous information and random effects

The regression framework allows to extend the model by including exogenous information. We consider again model (13) and parametrize the expectation through


with and again being the subject-specific sender- and receiver-effects. Furthermore, represents a -dimensional covariate vector and is the corresponding parameter vector. We can use the augmented Lagrangian approach from above to estimate the dimensional parameter vector . It is important to note here that only dyadic covariates have the potential to increase the predictive performance of the approach. If we only include subject-specific (monadic) information the expectation can be multiplicatively decomposed and the model collapses back to the IPFP model (14). Thus is easily seen through

Nevertheless, the inclusion of subject-specific information may be valuable if it is the goal to forecast future networks based on new covariate information. This holds in particular in dynamic networks. We give an example for predictions based on lagged covariates in the next section.

We can also easily add additional structure to model (17

) and assume a distributional form for some or all coefficients. A simple extension arises if we assume random effects. This occurs by the inclusion of normally distributed sender- and receiver-effects:


where we take as the vector of parameters that determines the covariance matrix of the random effects. The latter could be parametrized for example with such that


where we assume separate variance components for the sender- and the receiver effects, respectively. In order to fit the model, we follow a Laplace approximation estimation strategy similar to

Breslow and Clayton (1993). Details are given in the Appendix A.

4 Application

4.1 Data description

AT Austria ES Spain JP Japan
AU Australia FI Finland KR South Korea
BE Belgium FR France NL Netherlands
CA Canada GB United Kingdom SE Sweden
CH Switzerland GR Greece TR Turkey
CL Chile IE Ireland TW Taiwan
DE Germany IT Italy US United States of America
Table 1: Countries included in the analysis

The dataset under study is provided by the Bank for International Settlements (BIS) and freely available from their homepage. In general, the locational banking statistics (LBS) provide information about international banking activity by aggregating the financial activities (in million USD) to the country level. Within each country, the LBS accounts for outstanding claims (conceptualized as a valued dyad that consists of all claims banks from country to banks of country ) and liabilities of internationally active banks located in reporting countries (conceptualized as the reverse direction ). We have selected the 21 most important countries (see Table 1) for the time period from January 2005 to December 2017 as a quarterly series for the subsequent analysis. In Figure 1 the density of the network (number of existing edges relative to the number of possible edges) is shown on the left, the share of the zero-valued marginals in the middle and the development of the aggregated exposures on the right. Especially in the first years some marginals of the financial networks are zero and the corresponding matrix entries are therefore not included in the estimation problem. Correspondingly, it can be seen that most countries do have some claims and liabilities to other countries but especially in the beginning, many dyads are zero valued.

Figure 1: Density (left), share of zero-valued marginals (middle) and aggregated volume in million USD (right) of the network as quarterly time series.

Since it is plausible that financial interactions are related to the economic size of a country, we consider the annual Gross Domestic Product (GDP, in current USD Billions) as covariate. The data is provided for the years 2005-2017 by the International Monetary Fund on their homepage. Furthermore, there might be relationship between trade in commercial goods and financial transfers and we use data on dyadic trade flows (in current USD) between states as additional covariate. The data is available annually for the years 2005 to 2014 by the Correlates of War Project online. We do not have available information on trade for the years 2015, 2016 and 2017 and we therefore extrapolate the previous values using an autoregressive regression model. Apparently, by doing so we have covariate information which is subject to uncertainty. We ignore this issue subsequently. In order to have an example for uninformative dyadic information, we use time-invariant data on the dyadic distance in kilometres between the capital cities of the countries under study (provided by Gleditsch, 2013). Finally, in some matrix reconstruction problems, the matrix entries of previous time points become known after some time. Typically, lagged values are strongly correlated with the actual ones. We therefore also consider the matrix entries, lagged by one quarter as covariates. See Table 2 for an overview of the variables, together with the overall correlation of the actual claims and the respective covariate. In the subsequent analysis we include all covariates in logarithmic form.

Variable Description Type Correlation
Claim from country to country dyad specific
Gross Domestic Product of country node specific
Gross Domestic Product of country node specific
Bilateral trade flows of commercial goods dyad specific
Distance between the capital cities of countries and dyad specific
Lagged claim from country to country dyad specific
Table 2: Covariates used for the regression-based network-reconstruction

4.2 Model performance

Method Covariates Rand. eff. Model overall overall average SE average SE
1 IPFP - - (11,14)
2 Regression - (14,19)
3 Regression , , - (17) 3 300.242 56.794 63.466 7.802
4 Regression , , (17, 19 9.222 1.052
5 Regression , , - (17)
6 Regression , , (17, 19)
7 Regression , , - (17) 2 280.235 41.805 43.851 5.591 1.549
8 Regression , , (17, 19) 11.715
Table 3: Comparison of different regression models with the BIS Dataset. All values scaled by and lowest values in bold.

We evaluate the proposed models in terms of their and errors. The corresponding results are provided in Table 3. As a baseline specification, all models contain sender- and receiver effects. In the first row, we provide the maximum entropy model (14) that coincides with the IPFP solution (11). The second row shows model (14) together with the random effects structure (19). In the third row, we provide the errors for model (17) where we included the covariates logarithmic GDP (, ) as well as the logarithmic trade data (). In row four, we use the same model as in row three but additionally added the random effects structure from (19). In rows five and six, the same models as in rows three and four are used but with logarithmic distance () instead of trade as dyadic explanatory variable. In the last two rows we consider models with lagged claims () with and without random effects. This comparison might be somewhat unfair because of the strong correlation and because it is not clear whether it can safely be assumed that such data is always available. Therefore, we have separated this specification from the other models.

In the first four columns the different specifications together with the related equations are provided. Columns five and six show the aggregated errors over all 52 quarters and the last four columns show the errors averaged over all years together with their corresponding standard errors.

Figure 2: Quarterly time series of the estimated variances (models (14) and (19)) of the sender-effect () on the left, the receiver effect () in the middle and the correlation between the sender- and the receiver effect () on the right.

It can be seen that the first two models provide the same predictions and the inclusion of the random effects has no impact other than giving estimates for the variance of the sender- and receiver-effects as well as their correlation, shown in Figure 2. It becomes visible that the variation of the receiver-effect is much higher than the variation of the sender-effect which is almost constant. The correlation between the sender and the receiver effect is consistently positive and increases strongly within the first years.

Furthermore, Table 3 shows that in the four models that include exogenous information (rows three to six) the extension towards the random effects structure has an impact on the predictive quality. It decreases in the model that includes the variable and increases in the one that includes . Nevertheless, the models with and without random effects are rather close to each other and in fact they are statistically indistinguishable with respect to their differences. While the model with the covariate performs even worse than the IPFP solution, the model that includes but includes no random effects (row three) gives superior predictions relative to all other models in the upper part of Table 3. However, the two models that use the information on the lagged values give by far the best predictions. We nevertheless continue with the best model from the upper part of Table 3 lagged data are not necessarily available.

Figure 3: Estimated coefficients of model (17) with , and as covariates. Time series of sender- (, left) and receiver-effects (, right) in the first row. Time series of estimated coefficients on exogenous covariates () in the second row.

The corresponding fitted values are provided as time series in Figure 4 and in Figure 3 we provide the estimates for the coefficients of the model. In the first row, the estimated sender- (, left) and receiver-effects (, right) are shown as a time series. In the second row of Figure 3 the estimates for the coefficients on the exogenous covariates () can be seen. The estimated coefficients provide the intuitive result that the claims from country to country increase with and and the trade volume between them (). It is reassuring that the ordering of the average height of the coefficients approximately matches with the order of the correlations reported in Table 2. Note however, that the size of the coefficients is to be interpreted with care because of the limited information available on the unknown claims.

Figure 4: Fitted exposures, shown as quarterly series (52 time points) in million USD. All estimated with model (17) with , and as covariates. Four dyads with the overall highest values are highlighted.

We also provide prediction intervals in Figure 5, based on the share of real values located in the interval . Here, and are the and -quantiles derived from the bootstrap distribution (bootstrap sample size ). On the left, we illustrate the real values against the predicted ones together with grey 95% prediction intervals for the most recent network. Observations that do not fall within the prediction interval are indicated by red circles. Because of the quadratic mean-variance relation of the exponential distribution it is much easier to capture high values within the prediction intervals than low ones. A circumstance that materializes in the fact that exclusively small values are outside the prediction intervals. The share of real values within the prediction intervals against time is shown on the right hand side of Figure 5. We cover on average of all true values with our prediction intervals over all time periods and regard the bootstrap approach therefore as satisfying.

Figure 5: 95% Prediction intervals for the estimated means of model (17) with , and as covariates, for the most recent network on the left. Prediction intervals in grey. Real values on the horizontal axis and estimated means on the vertical axis. Bisecting line in solid black. Predictions (not) covered by prediction intervals in black (red) circles. Share of real values within prediction intervals as a time series on the right.

4.3 Comparison to alternative routines

Gravity model: A standard solution to the problem is the gravity model (e.g. Sheldon and Maurer, 1998). In essence it represents a multiplicative independence model


The model is simple, easy to implement and very intuitive. In situations where diagonal elements are not restricted to be zero it even coincides with the maximum entropy solution.
Tomogravity model: An extension of the gravity model is given by the tomogravity approach by Zhang et al. (2003)

. The model was initially designed to estimate point-to-point traffic volumes from dyadic loads and builds on minimizing the loss-function


subject to the non-negativity constraint. Here, the gravity model (20) serves as a null model in the penalization term and the strength of penalization is given by . The approach is implemented in the R package tomogravity (see Blocker et al., 2014). Zhang et al. (2003) show in a simulation study, that is a reasonable choice if no training set is available. In our competitive comparison we optimize the tuning parameter in order to minimize the overall error with grid search and find to be an optimal value.
Non-negative LASSO: The LASSO (Tibshirani, 1996) was already applied to predict flows in bike sharing networks by Chen et al. (2017) and works best with sparse networks. Using this approach, we minimize the loss function


where we can drop the absolute value in the penalization term because of the non-negativity constraint (see Wu et al., 2014 for the non-negative LASSO). In order to use the approach in the competitive comparison, we optimize the penalty parameter on a grid for the minimum error and use . The models are estimated using the R package glmnet by Friedman et al. (2009).
Ecological regression: In Ecological Inference (see e.g. Klima et al., 2016, King, 2013), it is often assumed that the observations at hand are independent realizations of a linear model, parametrized by time-constant transition-shares . Define the stacked column sums in by and the stacked row sums in by . Then, the model can be represented as


where the matrix contains the parameters . In the give form, the problem is not identified for one time period and it must be assumed that multiple time-points can be interpreted as independent realizations. Additionally, the model is not symmetric, implying that the solution to equation (23) does not coincide with the solution to


Since, the estimated transition shares are not guaranteed to be non-negative and sum up to one they must be post-processed to fulfil this conditions. Both models are fitted via least-squares in R.
Hierarchical Bayesian models: Gandy and Veraart (2017)

propose to use simulation-based methods. In their hierarchical models the first step consists of estimating the link probabilities and given that there is a link, the weight is sampled from an exponential distribution:


In order to estimate the link probabilities , knowledge of the density or a desired target density is needed. In their basic model it is proposed to use an Erdös-Rényi model with consistent with the target density. In an extension of the model, inspired by Graphon models, Gandy and Veraart (2018) propose a so called empirical fitness model. Here the link probability is determined by the logistic function


with being some constant that is estimated for consistency with the target density. For the fitness variables , the authors propose to use an empirical Bayes approach, incorporating the information of the row and column sums as . An implementation of both models is given by the R package systemicrisk. In order to make the approach as competitive as possible we use for each quarter the real (but in principle unknown) density of the networks. Because the results of the method differ between each individual estimate, we average the estimates and evaluate the combined dataset.

4.4 Competitive comparison

Method Model overall overall average SE average SE
1 Regression (, , ) (17) 3 300.242 56.794 63.466 9.246 7.802 1.085
2 Gravity model (20)
3 Tomogravity Model (21)
4 Non-negative LASSO (22)
5 Ecological Regression, columns (23)
6 Ecological Regression, rows (24)
7 Hierarchical, Erdös-Rényi (25)
8 Hierarchical, Fitness (26)
Table 4: Comparison of different methods with the BIS Dataset. All values scaled by and lowest values in bold.

Again we compare the different algorithms in terms of their and errors in Table 4. In the first row of Table 4 we show the restricted maximum likelihood model with the best predictions from Table 3 and in the following rows, the models introduced in Section 4.3 above are shown. The results can be separated roughly into three blocks. The models that fundamentally build on some kind of Least Squares criterion without referring in some way to the gravity model or the maximum entropy solution (ecological regression, and non-negative LASSO in rows four, five and six) have the highest values in terms of their and errors. Somewhat better are the Hierarchical Bayesian Models (rows seven and eight) that can be considered as the second block. However, although they provide better predictions than the models in the first block, we used the real density of the network in order to calibrate them which gives them an unrealistic advantage. The third group is given by the gravity and tomogravity model (rows two and three). Those are statistically indistinguishable and provide considerably better results than the models from the former blocks. Nevertheless, the regression model that uses exogenous information on (first row) yields the best predictions in this comparison.

4.5 Performance of the estimator

We hope to see improvements in the predictions if we include informative exogenous variables in the model. Informative means in this context, that variation in is able to explain variation in the unknown . Apparently, including information with a low association to simply adds noise into the estimation procedure. In this case we expect inferior predictions as compared to the IPFP solution. We illustrate the properties of the estimator using a simulation study with the following data generating process


Since the association between and the unknown is crucial, we vary the parameter from to and denote with the mean based on . For each parameter we re-run the data generating process (27) times and calculate for the -th simulation the IPFP solution and the restricted maximum likelihood solution . Based on that, we calculate in each simulation the ratio of the squared errors

This ratio is smaller than one if the IPFP estimates yield a lower mean squared error than the restricted maximum likelihood estimates and higher than one if the exogenous information improves the predictive quality in the terms of the mean squared error.

Figure 6: Median (solid black) and mean (solid grey) of the relative squared error for different levels of .

In Figure 6, we show the median (solid black) and the mean (solid grey) of for different values of as well as a horizontal line indicating the value one (dashed black) and a vertical line for (dashed black). It can be seen, that the mean and the median of are below one for values of that are roughly between and but increase strongly with higher absolute values of . Apparently, the distribution of

is skewed with a long tail since the mean is much higher than the median. With very high or low values of

, the median of the relative mean squared error becomes more volatile and partly decreases.

Figure 7: Boxplots for the normalized differences between the true expectation and the mean of the estimated ones for (top), (middle) and

(bottom). Regression-based model on the left, IPFP on the right. Outliers are excluded from the boxplot representation for better clarity.

Furthermore, we investigate how well the estimated expectations match the true ones. By construction, it holds that and consequently, the regression-based approach as well as IPFP assume that the sum of realized values equals the sum of the true expectations. Nevertheless, the moment condition does not imply that . In order to investigate potential bias, we draw again from

and fix to be the true expectation and draw and re-estimate again times from . Apparently, estimating the true expectations is a hard task as only sums of random variables with different expectations are available. Consequently, the variation of the mean estimates is rather high and we report boxplots of the normalized difference between the true value and mean estimate

and accordingly for . In Figure 7 we illustrate three different cases with (top), (middle) and (bottom). On the left hand side, boxplots for are shown for the regression-based model and on the right hand side for IPFP. The solid black line represents zero and the dashed black lines give 1.96. The results for the case on the top, match with the previous analysis illustrated in Figure 6 and show that IPFP identifies the true expectations somewhat better than the regression-based approach when the exogenous information is non-informative. In such a case, including adds noise in the estimation procedure, resulting in a greater variance around the true expectations. However, this changes strongly if is informative. Especially for on the bottom of Figure 7, some estimates obtained from IPFP are seriously biased because this procedure does not have the ability to account for the dyad-specific heterogeneity. The regression-based method, however does a reasonable job in recovering the unknown true expectations.

5 Discussion

In this paper we propose a method that allows for network reconstruction within a regression framework. This approach makes it easy to add and interpret exogenous information. It also allows to construct bootstrap prediction intervals to quantify the uncertainty of the estimates. Furthermore, the framework is flexible enough to deal with problems that involve partial information. For example if some elements of the network are known or if we have information on the binary network, then we can model the expected values of the matrix entries conditional this information, simply by changing the routing matrix and the E-Step.

However, we also want to list some shortfalls of the method. An obvious drawback of the method is its derivation from the maximum entropy principle that tries to allocate the matrix entries as even as possible and is therefore not very suitable for sparse networks as long as the sparseness cannot be inferred from the marginals. Furthermore, the estimated coefficients must be interpreted with care, as they are estimated based on a data situation with much less information than in usual regression settings. As a last but most important point, if the association between the exogenous explanatory variable(s) and the unknown matrix entries is low, the method is likely to deliver predictions that are worse than simple IPFP. It is therefore highly recommendable to use expert knowledge when selecting the exogenous dyadic covariates for regression-based network reconstruction.

Appendix A Estimation with random effects

In order to fit a model of the form


we follow a Laplace approximation estimation strategy similar to Breslow and Clayton (1993) and fix to some value and follow a Laplace approximation estimation strategy similar to Breslow and Clayton (1993), see also Tutz and Groll (2010). The moment condition implies a restriction for but not for and given some starting value , we can maximize the penalized log-likelihood (constant terms omitted)

subject to the moment condition . Therefore, the new optimization problem is given by


Define as the design matrix for the fixed effects and as the () random effects design matrix, this allows to write the generic mean as . Given that we have some estimate of , call it we can estimate the variance parameters with an approximation of the marginal restricted log-likelihood:


where , with and and consequently . The pseudo-observations are given by . Estimators can be obtained by iteratively optimizing firstly (29) and secondly (30) in each iteration until convergence.


Data Sets:

BIS bilateral exposures dataset (BIS_bilateral_exposures.csv), capital distances dataset (capdist.csv), dyadic trade flows (Dyadic_COW_4.0.csv), gdp values (gdp.csv)


A complete workflow ( that allows to reproduce the whole analysis with all relevant Figures and Tables. I.e. reading in the raw data in R, estimation routines for all models mentioned in the paper and functions to create relevant plots and tables.


  • Airoldi and Blocker (2013) Airoldi, E. M. and A. W. Blocker (2013). Estimating latent processes on a network from indirect measurements. Journal of the American Statistical Association 108(501), 149–164.
  • Anand et al. (2018) Anand, K., I. van Lelyveld, Ádám Banai, S. Friedrich, R. Garratt, G. Hałaj, J. Fique, I. Hansen, S. M. Jaramillo, H. Lee, J. L. Molina-Borboa, S. Nobili, S. Rajan, D. Salakhova, T. C. Silva, L. Silvestri, and S. R. S. de Souza (2018). The missing links: A global study on uncovering financial network structures from partial data. Journal of Financial Stability 35, 107 – 119.
  • Bacharach (1965) Bacharach, M. (1965). Estimating nonnegative matrices from marginal data. International Economic Review 6(3), 294–310.
  • Bacharach (1970) Bacharach, M. (1970). Biproportional matrices and input-output change. Cambridge: Cambridge University Press.
  • Bishop et al. (1975) Bishop, Y. M., P. W. Holland, and S. E. Fienberg (1975).

    Discrete multivariate analysis: theory and practice

    Cambridge: MIT Press.
  • Blocker et al. (2014) Blocker, A. W., P. Koullick, and E. Airoldi (2014). networkTomography: Tools for network tomography. R package version 0.3.
  • Breslow and Clayton (1993) Breslow, N. E. and D. G. Clayton (1993).

    Approximate inference in generalized linear mixed models.

    Journal of the American Statistical Association 88(421), 9–25.
  • Brown (1986) Brown, L. D. (1986). The dual to the maximum likelihood estimator. In G. S. Shanti (Ed.), Fundamentals of statistical exponential families with applications in statistical decision theory, Chapter 6, pp. 174–207. Hayward: Institute of Mathematical Statistics.
  • Chen et al. (2017) Chen, L., X. Ma, G. Pan, J. Jakubowicz, et al. (2017). Understanding bike trip patterns leveraging bike sharing system open data. Frontiers of computer science 11(1), 38–48.
  • Coates et al. (2002) Coates, A., A. O. Hero III, R. Nowak, and B. Yu (2002). Internet tomography. IEEE Signal processing magazine 19(3), 47–65.
  • Dahmen and Reusken (2006) Dahmen, W. and A. Reusken (2006). Numerik für Ingenieure und Naturwissenschaftler. Berlin, Heidelberg: Springer-Verlag.
  • Deming and Stephan (1940) Deming, W. E. and F. F. Stephan (1940). On a least squares adjustment of a sampled frequency table when the expected marginal totals are known. The Annals of Mathematical Statistics 11(4), 427–444.
  • Dudík et al. (2007) Dudík, M., S. J. Phillips, and R. E. Schapire (2007). Maximum entropy density estimation with generalized regularization and an application to species distribution modeling. Journal of Machine Learning Research 8(6), 1217–1260.
  • Dym and Shames (2013) Dym, C. L. and I. H. Shames (2013). Introduction to the calculus of variations. In Solid Mechanics, Chapter 2, pp. 71–116. New York: Springer.
  • Efron and Tibshirani (1994) Efron, B. and R. J. Tibshirani (1994). An introduction to the bootstrap. Boca Raton: CRC press.
  • Elsinger et al. (2013) Elsinger, H., A. Lehar, and M. Summer (2013). Network models and systemic risk assessment. In J.-P. Fouque and J. A. Langsam (Eds.), Handbook on Systemic Risk, Chapter IV, pp. 287–305. Cambridge: Cambridge University Press.
  • Fienberg et al. (1970) Fienberg, S. E. et al. (1970).

    An iterative procedure for estimation in contingency tables.

    The Annals of Mathematical Statistics 41(3), 907–917.
  • Friedman et al. (2009) Friedman, J., T. Hastie, and R. Tibshirani (2009). glmnet: Lasso and elastic-net regularized generalized linear models. R package version 2.0-1.6.
  • Gandy and Veraart (2017) Gandy, A. and L. A. Veraart (2017). A bayesian methodology for systemic risk assessment in financial networks. Management Science 63(12), 4428–4446.
  • Gandy and Veraart (2018) Gandy, A. and L. A. M. Veraart (2018). Adjustable network reconstruction with applications to CDS exposures. Journal of Multivariate Analysis, in press.
  • Gleditsch (2013) Gleditsch, K. S. (2013). Distance between capital cities. Accessed: 2017-04-07.
  • Golan and Judge (1996) Golan, A. and G. Judge (1996). Recovering information in the case of underdetermined problems and incomplete economic data. Journal of statistical planning and inference 49(1), 127–136.
  • Haberman (1978) Haberman, S. J. (1978). Analysis of qualitative data 1: Introductory topics. New York: Academic Press.
  • Haberman (1979) Haberman, S. J. (1979). Analysis of qualitative data 2: New Developments. New York: Academic Press.
  • Hazelton (2010) Hazelton, M. L. (2010). Statistical inference for transit system origin-destination matrices. Technometrics 52(2), 221–230.
  • Hestenes (1969) Hestenes, M. R. (1969). Multiplier and gradient methods. Journal of optimization theory and applications 4(5), 303–320.
  • Johnson (2014) Johnson, S. G. (2014). The NLopt nonlinear-optimization package. R package version 1.2.1.
  • King (2013) King, G. (2013). A solution to the ecological inference problem: Reconstructing individual behavior from aggregate data. Princeton: Princeton University Press.
  • Klima et al. (2016) Klima, A., P. W. Thurner, C. Molnar, T. Schlesinger, and H. Küchenhoff (2016). Estimation of voter transitions based on ecological inference: An empirical assessment of different approaches. AStA Advances in Statistical Analysis 100(2), 133–159.
  • Koller et al. (2009) Koller, D., N. Friedman, and F. Bach (2009). Probabilistic graphical models: principles and techniques. Cambridge: MIT Press.
  • Malvestuto (1989) Malvestuto, F. (1989).

    Computing the maximum-entropy extension of given discrete probability distributions.

    Computational Statistics & Data Analysis 8(3), 299–311.
  • Mastrandrea et al. (2014) Mastrandrea, R., T. Squartini, G. Fagiolo, and D. Garlaschelli (2014). Enhanced reconstruction of weighted networks from strengths and degrees. New Journal of Physics 16(4), 043022.
  • Meng and Rubin (1993) Meng, X.-L. and D. B. Rubin (1993). Maximum likelihood estimation via the ECM algorithm: A general framework. Biometrika 80(2), 267–278.
  • Miller and Blair (2009) Miller, R. E. and P. D. Blair (2009). Input-output analysis: foundations and extensions. Cambridge: Cambridge University Press.
  • Muñoz-Cobo et al. (2017) Muñoz-Cobo, J.-L., R. Mendizábal, A. Miquel, C. Berna, and A. Escrivá (2017).

    Use of the principles of maximum entropy and maximum relative entropy for the determination of uncertain parameter distributions in engineering applications.

    Entropy 19(9), 486.
  • Onuki (2013) Onuki, Y. (2013). Extension of the iterative proportional fitting procedure and its evaluation using agent-based models. In M. Tadahiko, T. Takao, and T. Shingo (Eds.), Agent-Based Approaches in Economic and Social Complex Systems VII, pp. 215–226. Berlin, Heidelberg: Springer.
  • Powell (1969) Powell, M. J. (1969). A method for nonlinear constraints in minimization problems. In R. Fletcher (Ed.), Optimization, pp. 283–298. New York: Academic Press.
  • Schneider and Zenios (1990) Schneider, M. H. and S. A. Zenios (1990). A comparative study of algorithms for matrix balancing. Operations Research 38(3), 439–455.
  • Sheldon and Maurer (1998) Sheldon, G. and M. Maurer (1998). Interbank lending and systemic risk: An empirical analysis for switzerland. Swiss Journal of Economics and Statistics 134(4), 685–704.
  • Squartini and Garlaschelli (2011) Squartini, T. and D. Garlaschelli (2011). Analytical maximum-likelihood method to detect patterns in real networks. New Journal of Physics 13(8), 083001.
  • Tibshirani (1996) Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 58(1), 267–288.
  • Tutz and Groll (2010) Tutz, G. and A. Groll (2010). Generalized linear mixed models based on boosting. In Statistical Modelling and Regression Structures, pp. 197–215. Berlin, Heidelberg: Springer.
  • Upper (2011) Upper, C. (2011). Simulation methods to assess the danger of contagion in interbank markets. Journal of Financial Stability 7(3), 111–125.
  • Vardi (1996) Vardi, Y. (1996). Network tomography: Estimating source-destination traffic intensities from link data. Journal of the American Statistical Association 91(433), 365–377.
  • Wu et al. (2014) Wu, L., Y. Yang, and H. Liu (2014). Nonnegative-lasso and application in index tracking. Computational Statistics & Data Analysis 70, 116–126.
  • Zhang et al. (2003) Zhang, Y., M. Roughan, C. Lund, and D. Donoho (2003). An information-theoretic approach to traffic matrix estimation. In Proceedings of the 2003 conference on Applications, technologies, architectures, and protocols for computer communications, pp. 301–312. Association for Computing Machinery.
  • Zhou et al. (2016) Zhou, H., L. Tan, Q. Zeng, and C. Wu (2016).

    Traffic matrix estimation: A neural network approach with extended input and expectation maximization iteration.

    Journal of Network and Computer Applications 60, 220 – 232.