Asymmetric Relatedness from Partial Correlation

Relatedness is a key concept in economic complexity, since the assessment of the similarity between industrial sectors enables policymakers to design optimal development strategies. However, among the different ways to quantify relatedness, a measure that takes explicitly into account the time correlation structure of exports is still lacking. In this paper, we introduce an asymmetric definition of relatedness by using statistically significant partial correlations between the exports of economic sectors and we apply it to a recently introduced database that integrates the export of physical goods with the export of services. Our asymmetric relatedness is obtained by generalising a recently introduced correlation-filtering algorithm, the partial correlation planar graph, in order to allow its application on multi-sample and multi-variate datasets, and in particular, bipartite temporal networks. The result is a network of economic activities whose links represent the respective influence in terms of temporal correlations; we also compute the statistical confidence of the edges in the network via an adapted bootstrapping procedure. We find that the underlying influence structure of the system leads to the formation of intuitively-related clusters of economic sectors in the network, and to a relatively strong assortative mixing of sectors according to their complexity. Moreover, hub nodes tend to form more robust connections than those in the periphery.


Introduction
In the past few years, the use of bipartite networks for the representation of realworld complex systems has become widespread in a variety of fields and applications. These networks are usually constructed using multi-sample, multi-variate structured data used to model complex systems such as biological networks (enzymes and reactions [1], genes and diseases [2], plants and pollinators [3]), movies and actors [4,5], authors and papers [5,6], board of directors members and companies [7,8], companies and technologies they patent [9], members of peer-to-peer networks and data provided [10], international NGO branches and cities hosting them [11], supreme court judges and their votes [12], and legislators and bills they sponsor [13].
A prominent example is the bipartite network formed by countries and the products they export. This type of data has been used extensively in the field of economic complexity (EC) [14,15] to assess various quantities of interest for the modelling of the economic development of countries. The first one is the competitiveness of countries and the sophistication of products [16][17][18][19][20], and the relatedness between products, countries, or between countries and products [21,22]. With respect to the datasets implemented in the literature up to now, the dataset we use in this paper adds the inclusion of services to the set of tangible products traditionally considered in the EC literature [23,24]. by countries, averaging these, and applying the existing PCPG procedure. Following similar principles, we also adapt an existing bootstrapping procedure (see [34]) in order to determine the statistical reliability of the links present in the resulting network.
The contribution of this paper is many fold. Firstly, the biPCPG framework opens the possibility of the application of the PCPG algorithm to a wide variety of datasets with a multi-sample and multi-variate structure, including, but not limiting to, the ones usually analysed using the EC framework. Furthermore, the data-processing methodology introduced here could be utilised to apply other correlation-filtering algorithms for network generation (e.g., [31,33]).
Secondly, this paper introduces a network which describes the asymmetric relatedness among physical products (manufacturing) and services. This is an addition with respect to the networks usually present in the literature, such as the product space [21] and product taxonomy network [22], which are constituted only by products.
Thirdly, this paper introduces an adapted bootstrapping procedure to asses the reliability of the edges present in a network generated from multi-sample multi-variate datasets. Similarly to the network-generating framework, this bootstrapping procedure can be utilised to asses the reliability of edges in networks generated using alternative correlation-filtering methods with datasets with this structure.
Fourthly, in order to assess the information content of the biPCPG network we calculate two assortativity measures and run a community detection procedure, finding that meaningful clusters and connections emerge, as well as a relevant complexity-related assortativity. In summary, the biPCPG analysis unveils the average influence between industrial and service sectors, efficiently encapsulating the information about the correlation structure of the system.
Finally, we provide a Python package named "biPCPG" [35] with its documentation hosted in [36]. The 0.1.0 version of this package was used to perform all the calculations done in this paper, including the data-handling, biPCPG network generation, bootstrapping procedure and calculations done on the biPCPG network. It is worth noting that the package has a modular structure such that the data-handling and the generation of the biPCPG network are computed independently of each other. This allows the user to, for example, utilise the data-handling module to prepare a multi-sample multi-variate dataset for an alternative correlation-filtering method, or to implement the PCPG algorithm on a dataset of her choice, without the need for the dataset to have a multi-sample multi-variate structure. To the best of our knowledge, the PCPG module in the biPCPG package is the first publicly available Python implementation of the PCPG algorithm.
The rest of this paper is organised as follows. In Section 2, we describe the dataset used in this investigation and the cleaning procedure performed on it. In Section 3, we describe the set of methods to generate the biPCPG network and comment on the resulting network. In the result sections we describe the assortativity calculations and community detection procedure done on the biPCPG network and show the results obtained. Section 5 concludes.

Data Description and Preprocessing
The dataset used in this research project is an integration of the United Nations Commodity Trade Statistics Database (UN-COMTRADE-https://comtrade.un.org, accessed on 13 February 2019) and the International Monetary Fund's Balance of Payments data (BPM6) [37], relative to physical goods and service exports respectively. This integrated dataset was introduced in a World Bank working paper [23]. The UN-COMTRADE data consists of the amount of exports from each country per category of products (in USD). The categorisation of products is given by the World Customs Organization's (WCO) Harmonized System 2007 edition (HS2007) [38], which classifies products by using a hierarchical six-digit code depending on the category of the product. The IMF BPM6 dataset consists of the amount (in USD) of services provided abroad by each country and is collected according to the 6th edition of its manual, provided by the International Monetary Fund (IMF). Henceforth, we will globally refer to the collection of products in COMTRADE and services in BPM6 as sectors.
The hierarchical structure of the HS classification allows for an aggregation from the most granular six-digit level, consisting in about 5000 different products, into a coarser twodigit level. A further aggregation of a few small (in terms of export quantities) two-digit sectors into a single two-digit sector was also performed in this dataset, leaving a total of 78, roughly homogeneous aggregated product sectors at the two-digit level. From the BPM6 part of the dataset, there are a further 22 service sectors at a comparable level of aggregation.
The aggregated dataset used in our study is therefore comprised of 78 + 22 = 100 sectors of products and services, these are listed in Table A1 in Appendix D. The data span a total of 22 years, from 1995 to 2016. As there are missing data points in some years for several countries, we apply a sanitation procedure where only countries with complete data for all sectors throughout the 22 years are kept. This reduces the dataset to from 129 countries to 99 countries. The analysed dataset has a total 99 × 100 = 9900 time series of length 22, with no missing values, representing the amount of product exports or service provisions in USD for each country.
In order to perform specific calculations (see Section 4.2), the 100 sectors in the dataset must be aggregated one level further. The product sectors can be further aggregated using what the WCO refers to as sections. The WCO provides a total of 21 sections which are available at [38]. In this case, services sectors can be aggregated into a single "section". Thus, in our aggregated dataset we have a total of 22 sections of sectors-21 product sections arising from the HS2007 classification, and one additional section containing the service sectors from the BPM6 dataset.

Revealed Comparative Advantage Matrices
The raw data used to construct in this paper are the amount of exports E y c,p (in USD) of a sector p (product or service) by a country c in year y. We compute the Revealed Comparative Advantage (RCA) [39] as where P and C are the sets of unique sectors and unique countries in the dataset discussed above.
The use RCA is ubiquitous in the EC literature, because removes trivial dependencies from the sectors' and countries' size. When the RCA y c,p is above 1, the country is said to have a revealed comparative advantage in exporting a given sector in that year. Conversely, when RCA y c,p is below 1 the country can be thought of as not being very competitive in that particular sector. Finally, when RCA y c,p is equal to 1 the country has the expected (average) share of the world's exports in the given sector and year. Therefore, the dataset on which we perform the following calculations consists of time series RCA c,p = (RCA y c,p : y ∈ Y) for 99 countries and 100 sectors, where Y is the index set of years [1995,2016]. The data is then shaped into a set of 22 matrices RCA y , one for each year, where each row represents a country, each column represents a sector and each entry is the corresponding RCA y c,p value.

Methodology Description
Before discussing the detailed implementation of the biPCPG methodology, here we provide a summarised description of our procedure; a visual representation can be found in Figure 1.
Given the multi-sample nature of the dataset analysed, a series of data-preprocessing steps are needed before the application of PCPG. The PCPG algorithm takes a single correlation matrix as an input and outputs a network (see Section 3.5). In order to obtain our biPCPG network, along with reliability values for its edges from a multi-sample dataset, we need two main procedures, a "Network generating procedure" and a "Bootstrapping procedure".
The "Network generating procedure" is shown in the black box in Figure 1 and deals with the data handling necessary to obtain a PCPG network from a dataset with a multisample structure. In our case, we are interested in obtaining a biPCPG network where nodes are sectors, therefore the input matrix should describe the correlations between sectors.
To find this input correlation matrix, the initial step is to shape the dataset such that, for each country, we have a matrix where the columns are the relevant time series of each sector. We then compute a correlation matrix for each of these time series matrices. Finally, we average these correlation matrices over countries to obtain an average correlation matrix which serves as the input to the PCPG algorithm, i.e., the last step in the biPCPG framework. The output of the biPCPG algorithm is the network we refer to as G, as well as the weights of the edges in contains, i.e., the average influence between sectors. The "Bootstrapping procedure" of our framework, shown in the grey box in Figure 1, deals with the bootstrapping procedure necessary to asses the reliability of the edges in the biPCPG network obtained. This starts from the country time series matrices, which are bootstrapped R times, obtaining a "batch" of replicates each time. Each of these batches contains C matrices, one for each country, where the rows have been drawn coherently from their corresponding original country matrices. This is done in order to randomise the time dimension while preserving the correlation structure across countries (see Section 3.6). We then replicate the "Network generating procedure" described above by treating each batch of replicates as a new dataset of country time series matrices and follow the steps to obtain a replicate biPCPG network. This means that, for each batch, we calculate a correlation matrix for every time series matrix, we then average across these correlation matrices and use the average correlation matrix as an input to the PCPG algorithm. Repeating this procedure for all R batches we obtain R replicate networks. We find the fraction of times each edge in G appears in the replicate networks, which is a measure of the reliability of the edge.

Partial Correlations and Average Influence: Definitions
As described in the original PCPG paper (see [30]), the starting point of our analysis is the partial correlation, which measures the effect that a random variable Z has on the correlation between two other random variables, X and Y. The partial correlation ρ(X, Y : Z) is defined in terms of the Pearson correlations ρ(·, ·) between the three variables, formally A small value of ρ(X, Y : Z) may be ambiguous, as this could be due to the correlations among the three variables being small; or due to variable Z having a strong effect on the correlation between X and Y, which is generally the interesting case. In order to discriminate between these two cases the correlation influence or influence of variable Z on the pair of elements X and Y is used. This is defined as We define the average influence of variable Z on the correlations between X and all other variables in the system as follows: We anticipate that the average influence will be the input of the network building algorithm also described in [30]. Note that, potentially, there could be certain values of measured correlations ρ(X, Y), ρ(X, Z) and ρ(Y, Z) that lead to a measured partial correlation ρ(X, Y : Z), to be out of its defined range [−1, 1]. In our analysis, this occurred in 0.02% of the partial correlations computed. In these cases, partial correlations were set to be undefined (NaN in programming terms) which in turn makes the influence values based on these partial correlations also undefined. Similarly to the undefined correlation values described above, these undefined influences are not included in calculation of average influence d(X : Z).
Some of the values obtained for ρ(X, Y), ρ(X, Y : Z), d(X, Y : Z) and d(X : Z) in our dataset and their interpretation are discussed in Section 3.4. An important point is that, in general, d(X : Z) = d(Z : X): the influence is asymmetric, and the largest among these two quantities indicates the main direction of influence between X and Z. For example, in our dataset when X = Glass and Z = Furniture, the average influence of Furniture on Glass d(X : Z) = 0.03 while the corresponding reverse average influence of Glass on Furniture d(Z : X) = 0.29, suggesting that the direction of influence is from Glass to Furniture and not vice-versa. This, however, is an example of a clear-cut case, where difference between the two average influence values is not small. In general, these differences tend to be much smaller. This can be an effect of the complex relationship and mutual interaction between the economic sectors, or a consequence of the noise present in the data. This makes a bootstrapping procedure necessary in order to asses the statistical confidence in the overall direction of influence, as well as the average influence values themselves. We will discuss the bootstrapping procedure in Section 3.6.

Average Correlation Matrix
The input to the PCPG algorithm is a correlation matrix [30]. In our procedure, to allow its use on our multi-sample dataset, this correlation matrix is replaced by an average correlation matrix over countries. In order to obtain this average correlation matrix, we reshape the 22 RCA y matrices into a total of C = 99 matrices, one for each country, each consisting of T = 22 rows and P = 100 columns. We denote these TS c , c ∈ 1, . . . , C.
In this way, the columns of each matrix TS c are the RCA c,p time series of the given country c, where each column represents a sector p in the dataset.
In order to obtain the input matrix to the PCPG algorithm, we first find C correlation matrices denoted K c , c ∈ 1, . . . , C from the pair correlations between the columns of each matrix TS c . Thus the entries of the country correlation matrix K c are given by where ρ is the Pearson correlation, the subscript * , p denotes the column p of the matrix and RCA c,p is the RCA time series for country c and sector p.
For each correlation value we obtain p-value via a two-sided T-test procedure [40]. Given we are performing multiple tests, we apply a False Discovery Rate (FDR) correction to obtain adjusted p-values via the Benjamini-Hochberg (BH) procedure [41]. We choose the BH procedure since it ultimately allows the inclusion of more information in the biPCPG network than a more restrictive correction procedure such as the Bonferroni correction [42]. Note that the FDR correction has been extensively used in the literature for the statistical validation of networks and, in particular, it has been previously used to validate networks representing bipartite complex systems [43].
We reject non-statistically significant correlation samples when the adjusted p-value is above a critical value of 0.01. In these cases, the corresponding entries to the K c matrix are marked as undefined. The same procedure for obtaining country correlation matrices was also performed without the FDR correction for the 0.01 and alternative critical values. This produced networks which have the same main features as the network presented below, including the main hub nodes, clusters of sectors and communities detected.
Once the country correlation matrices K c are found, we then compute the element-wise mean of these matrices, obtaining the average correlation matrixK with entries where row and column indices p and p denote economic sectors. Any undefined correlation is discarded during the averaging process. Note that, using this notation, the correlations ρ(·, ·) mentioned in Section 3.2, are replaced by the average correlationsK p,p described here. This leads to an equivalent expression for the partial correlation

Partial Correlation and Average Influence: Empirical Analysis
In order to clarify the meaning of the intermediate quantities that are used to build the biPCPG network, we devote this subsection to the discussion of some empirical features.
Bearing in mind how the influence of a variable on the correlation of two other variables is defined (see Equation (3)), we explore four examples of the results obtained from these computations. Note that, in the description below, the variables X, Y and Z used in the definition of Equation (3), are replaced by sectors of our system. Thus, the partial correlation column in Table 1 describes the average correlation,K p,p , between sectors p and p accounting for the effect of a third sector p , and similarly for the influence column. We therefore denote these quantities ρ(p, p : p ) and d(p, p : p ), respectively. Example 1 shown in Table 1 is an example of the case described in Section 3.2, which shows a very small partial correlation due to all correlations among the three variables being small. By definition, this makes the resultant influence value is small, which reduces the average influence of the sector "Other textile" on the sector "Cereals", making the appearance of this edge in the network less probable.
Example 2 also shows a case where the partial correlation between p and p , accounting for the effect of p , is small. However, contrary to the case in Example 1, this is due to p strongly affecting the correlation between p and p , i.e., ρ(p, p ) ∼ ρ(p, p )ρ(p , p ). Therefore, the resulting influence is relatively high, which increases the probability of an edge from "Cultural" to "Audiovisual" being present in the biPCPG network. In addition, note that the probability of an edge from "Cultural" to "Audiovisual" also increases under these results, due to the symmetry between the p and p variables.
In Example 3, we have a case where the correlation between p and p is relatively strong and variable p has a small effect on it. This is due to the similar values of the correlation ρ(p, p ) and the partial correlation ρ(p, p : p ). Therefore, the resulting influence of "Knitted clothing" on the correlation between the "Pigments" and "Aluminium" sectors is close to zero.
Finally, Example 4 shows a seemingly counter-intuitive case where the correlation between p and p is small while their partial correlation given p is negative, yielding a high influence. A negative partial correlation occurs when the correlation between p and p is small but both p and p have a high correlation with p . In this case, the influence of "Plastics" can be interpreted as preventing the correlation ρ(p, p ) between "Vehicles" and "Earths and stone" from being lower, or being negative. It is important to note that the average influence values among sector pairs determine the structure of any PCPG network (see Section 3.5). Figure 2 displays a scatter plot that shows the correlation ρ(·, ·) and average influence d(·, ·) among all N(N − 1) = 9900 pairs of sectors in our biPCPG network. Note that this includes data points for both d(p : p ) and d(p : p) influences at the same horizontal coordinate as the correlation between p and p is symmetric.
This plot shows that the average influence between a pair of sectors is highly correlated with the correlation between the same pair of sectors, showing a very narrow 95% confidence interval (barely visible as it is only slightly wider than the fit line). See Appendix B for details on the calculation of the confidence and prediction intervals shown in Figure 2. This is not surprising given how the average influence is calculated; however, the relatively high coefficient of determination R 2 = 0.58 indicates that, generally, the partial correlation values obtained are relatively small. This may be due to there actually not being large influences between the sectors, or due to limitations of the dataset. For example, hidden influences between the sectors could potentially be detected in datasets with longer time series.
In Figure 2, we can observe that most of the correlations (around 80%) are positive. Around 10.7% of the pairs of sectors with positive correlations have an average influence below zero. This quantity is over an order of magnitude larger than its counterpart, the percentage of pairs of sectors with negative correlation but a positive average influence, which is around 0.47%.

Figure 2.
Plot showing correlation and average influence values among all 9900 pairs of sectors in the system. A line of best fit among the points is shown in red along with the coefficient of determination R 2 = 0.58, with the 95% confidence interval limits in light blue and the 95% prediction interval limits in dashed grey lines. Note the confidence interval is so narrow it is only visible at the edges of the red best fit line upon close inspection.

Network Construction
The construction algorithm of a PCPG network starts with a list of the N(N − 1) average influence values in decreasing order and an empty graph of N nodes and no edges, where N is the number of variables in the system. In our case, we have N = 100 economic sectors. We then cycle through the sorted list, starting with the largest average influence value found, e.g., d(p : p ), where p and p are a given pair of products. The edge p → p is included in the network if and only if the resulting network is still planar and the edge p → p has not been included already. We stop adding edges if adding the next edge in the list would break the planarity of the graph. This procedure ensures two things: (i) only the largest among d(p : p ) and d(p : p) will be included in the network, and (ii) the final network has 3(N − 2) edges. It is important to note that for a given input correlation matrix of size N × N the PCPG network will always have 3(N − 2) edges and that the identity of these edges solely depends on the correlation values in the input matrix.
The final result of this procedure is what we refer to as the biPCPG network, G. Naturally, we also obtain the average influence d associated to each edge in G, as well as the network's adjacency matrix A defined as

biPCPG Bootstrapping
To assess the reliability of the links in the biPCPG network, we adapt a bootstrapping procedure originally introduced in [34]. The aim is to obtain a bootstrap value for each link which is proportional to the reliability of the link.
We build R batches, where the matrices to be bootstrapped in each batch are the time series matrices of all countries TS c ∀ c ∈ 1, . . . , C. From each matrix TS c , a replicate time series matrix TS r c ∀ r ∈ 1, . . . , R is obtained, where R = 1000 is the total number of batches. An important feature of our procedure is how the null model, i.e., the replicate time series matrices, is generated. For each batch, the bootstrapping of the time series matrices is done coherently across countries. This means rows are drawn with repetition from each of the country matrices jointly-the same row indices are selected across the matrices. In addition, the new locations of the selected rows in their corresponding replicate matrices are exactly the same. This way, in the replicate time series matrices, TS r c , the time structure of the time series is destroyed while preserving the country-level correlations.
Take, for example, the first batch, r = 1. In order to obtain the first batch of replicate matrices TS 1 c ∀ c ∈ 1, . . . , C, we randomly select a sequence of T = 22 row indices, allowing repetitions. These row indices denote which rows from the original matrices TS c are included in the corresponding replicates TS 1 c in this batch, as well as their order. This way, any row of a replicate matrix in this first batch will contain data points corresponding to the same year as rows of the same index in all the other replicate matrices in the batch.
After all the replicate matrices are obtained for all countries and batches, we calculate a replicate correlation matrix K r c for each of them, rejecting non-statistically significant samples as described in Section 3.3. We then find the element-wise mean of the replicate correlation matrices in each batch r, obtaining R replicate average correlation matrices K r whereK Note that, similarly to the replicate time series matrices, in these replicate correlations matrices the time structure of the time series is destroyed while preserving the countrylevel correlations due to the way the bootstrapping has been performed. We then apply the PCPG algorithm described in Section 3.5 to each matrixK r , obtaining R replicate adjacency matrices, A r ∀ r ∈ 1, . . . , R.
To compute the bootstrap value, b p,p , for each link p → p , we evaluate the number of time the link appears in the replicate adjacency matrices A r , and normalise by the number of replicates R, formally Each bootstrap value is therefore some number in the interval [0-1] and is proportional to the reliability of the link.

Descriptive Analysis of the biPCPG Network
The network G resulting from the application of the biPCPG method to our dataset is shown in Figure 3. This network displays some interesting results with a few distinct hub nodes. The most noticeable of these nodes are "Plastics", "Pigments" and "Vegetables" nodes. Hub nodes in the network also tend to have high average influence on other nodes in the network, this being displayed by the width of the edges stemming out of them. The colour of the edge represents its bootstrap value. We note that the hub nodes are also the source of most of the darker edges in the network, i.e., the most reliable edges, especially the "Plastics" node, whose edges bootstrap values are very high. Figure 3. The biPCPG network. The widths of the edges are proportional to the average influence value, d(p, p ) they represent. The colours of the edges are proportional to their bootstrap value, b p,p . The darker the edge, the more reliable it is. Node colours represent the sector section each product and service belong to. Node sizes are proportional to out-degree. The node layout was found using the ForceAtlas2 algorithm [44].
The resulting network also displays distinct clusters of intuitively related economic sectors. For example, the most recognisable "food and plant" cluster can be found at the bottom-right of the network, surrounding the "Vegetables" hub node. At the top-left of the network, we can observe another distinct cluster containing several sectors related to chemicals or raw materials. Finally, on the top-right of the network, surrounding the "Plastics" and "Pigments" nodes, one can find a "macro-cluster" formed mostly by industrial and manufacturing sectors.
It is worth noting that, while most edges connect intuitively related sectors, the are several cases of less-intuitive connections spread around the network. This causes the inclusion of some of these seemingly unrelated sectors in some of the clusters mentioned above. This is partially due the original construction of the PCPG algorithm, which ensures a fixed number of edges to be included in the network. Therefore, edges representing small influences among sectors could be forced to be included in the network. In our case, the biPCPG network obtained contains around 5% of edges representing Average influence values of 0.05 or smaller.

Assortativity Analysis
As described in Section 2, the 100 sectors in our dataset can be grouped into 22 groups of sectors called sections. Furthermore, a key metric within the field of economic complexity is the complexity of a product or service, which measures the capabilities needed by a country to produce it (see Appendix A). In order to better understand the structure of this network, and by extension the information contained in it, one can then investigate its homophily or assortativity according to these characteristics. Roughly speaking, this is the tendency for nodes belonging to the same group to be connected to each other. In this paper, we make use of two different assortativity metrics which we describe below. The motivation behind this analysis is to assess if our framework generates a meaningful network which is able to synthesise information about the system.

Assortativity by Unordered Characteristics
This quantity is used to measure the assortativity between, for example, nodes with an associated qualitative characteristic such as, in our case, sector sections, s (see Section 2). The assortativity coefficient is defined as [45] where entries of the matrix F are the fractions of edges in the network that connect a vertex of section s to one of section s , and ||X|| is the sum of all elements of a matrix X [45]. Therefore the numerator is a quantity that measures the fraction of the edges in the network that connect vertices of the same type (i.e., within-section edges) minus the expected value of the same quantity in a network with the same community divisions but random connections between the vertices. The denominator is one minus the same expected value. This formula gives s s = 0 when there is no assortative mixing and s s = 1 when there is perfect assortative mixing. For a perfectly disassortative network, the value is in the range −1 ≤ s s < 0 (see [45] for its interpretation). We evaluate this metric for the section of sectors described in Section 2, denoting this by the subscript s.

Assortativity by Scalar Characteristics
A measure of assortativity for numeric quantities associated with nodes can also be defined [45]. First, note that the entries of the matrix F are the fraction of all edges in a network that connect nodes with associated scalar values q and q . Note that the values q and q are discrete-in our case these are the Complexity rank [17] of sectors-computed by taking average complexity value of each product (across the available years in our dataset) and ranking these averages from highest to smallest. The complexity of a product or service is a well-known quantity in the economic complexity literature that describes the capabilities needed by a country to produce it, see Appendix A for its definition. The numeric assortativity coefficient is defined as where a q = ∑ q F q,q , b q = ∑ q F q,q and σ a and σ b are the standard deviations of the distributions of a q and b q , respectively. The value of s q is in the range −1 ≤ s q ≤ 1 with s q = 1 indicating perfect assortativity and s q = −1 indicating perfect disassortativity. Typically, assortativity values in the range 0.3-0.7 are considered to indicate a significant community structure in social networks (higher values are rare) [46,47].

Assortativity Results
The results for the two assortativity metrics defined above are as follows: • assortativity by sector section = s s = 0.08 (0.15 without FDR correction); • assortativity by sector mean complexity rank = s q = 0.19 (0.31 without FDR correction).
These results indicate that the structure of the resulting biPCPG network encodes information efficiently. Firstly, the Assortativity by sector section, s s = 0.15, is positive, this means that sectors that belong to the same section (see Section 2) tend to be connected in the network, i.e., they influence each other. The section of each sector is reflected in Figure 3 by the colour of the node. The most evident clustering of sectors within the same section is found at the top of the plot where a highly connected cluster of service sectors is found.
Furthermore, the moderately high Assortativity by sector mean complexity rank, s q = 0.19, indicates that sectors around the same level of complexity tend to influence each other. This makes sense intuitively since, according to the economic complexity literature, these tend to be connected in other networks that describe the relationship among products (e.g., product space network, product taxonomy network [21,22]).

Community Detection on the biPCPG Network
We apply a well-known community detection algorithm for directed networks based on spectral optimisation [48]. The modularity, or quality function, to be maximised is where A is the adjacency matrix, k in p and k out p are the weighted in-degree and out-degree of node p, m is the total edge weight in the network, ν p is the community of node p and δ ν p , ν p = 1 if ν p = ν p and 0 otherwise. This method does not require any parameter choices relating to community size or number of communities; however, adaptations of this method that allow for these choices are available in the literature. It is worth pointing out that, for the analysis carried out in this paper, edge-weights are all set to 1. In Equation (13), this makes the weighted in-degree and out-degree simply the in-and out-degree as well as fixing m = 294, the total number of edges in the network.
Since there is no universal definition for communities in directed networks, we also apply the same community detection algorithm for the undirected version of the biPCPG network G und . In this case, the modularity to be maximised is given by where A und is the undirected adjacency matrix which defines the undirected network G und . This can be obtained from the adjacency matrix, A, which defines the directed biPCPG network G as follows This allows us to qualitatively assess if the structure of the biPCPG network is sufficient for reasonable communities to be detected, without the bias of the information contained in the average influence or bootstrap values associated to edges. We implement this algorithm via the leidenalg Python package (version 0.8.4) [49], an implementation of the leiden algorithm for modularity optimisation.
Note that optimising modularity is an NP-hard problem [50], and therefore heuristics have to be implemented for algorithms to be efficient. One of the steps in the leiden algorithm used here involves selecting a random community for a node to be added to. However, this randomness can be controlled via a seed to the random number generator. This makes the process deterministic such that the same communities are selected every time the algorithm is run on a given network using the same seed value. In our analysis, we tested several seed values finding that the detected communities varied only for a few nodes, with many seed values returning the exact same partitions. The results shown in Section 4.3 were found using 1 as the seed, as well as for many other seed values tested.
Furthermore, we compare the the communities obtained for the directed and undirected versions of the network for seed values 1, . . . , 1000 via the Adjusted Mutual Information [51]. Take, for example, our set of P of N sectors and consider two partitions of P, namely U = {U 1 , U 2 , . . . , U J } with J pairwise-disjoint clusters found by maximising Q und for the undirected version of the network, and V = {V 1 , V 2 , . . . , V D } with D pairwise-disjoint clusters found by maximising Q dir for the directed version of the network. The AMI between the two partitions is then defined as where MI(U, V) is the mutual information between two partitions, E{MI(U, V)} is the expected mutual information and H(U) and H(V) are the entropy values associated to partitions U and V respectively. The AMI equals 1 when two partitions are exactly the same and 0 when the MI between them equals its expected value and therefore serves as a similarity measure for the two partitions, for further details on its calculation see [51]. In Section 4.3, we give the result for the average AMI obtained for the 1000 seed values tested using the scikit-learn 0.23 Python package.

Community Detection Results
The community detection procedure described above yielded 5 distinct communities when applied on the undirected biPCPG network, G und , which we denote communities ν = 1, . . . , 5. These communities have 31, 22, 21, 13 and 13 sectors contained in each of them, respectively.
The detected communities in the network can be seen highlighted in Figure 4. When comparing with Figure 3, which shows the network highlighting the section of each sector, one can see that the detected communities partition the network into groups that contain intuitively related sectors. For example, communities 2, 3 and 5 contain mostly nodes related to industrial and chemical sectors, while community 1 captures the "food and plant" cluster described above as well as some service sectors. Finally, for community 4, it is slightly more difficult to find a common theme. However, it is worth noting that over half of the sectors it contains are service sectors. Figure 4. biPCPG network, G, resulting from the application of the PCPG algorithm on the mean correlation matrixK between sectors' RCA time series. Nodes are grouped by their community, ν, found by maximising modularity in the network. The node layout was found using the ForceAtlas2 algorithm [44].
The information structure these communities contain can be seen when sorting rows and columns of the average correlation matrixK and average influence matrix by community index as seen in Figures A2 and A3 in Appendix C. We can observe, for example, that brighter colours, meaning higher values, are generally found close the diagonal of the matrices (i.e., among sectors within the same community). This is especially noticeable for communities 1 and 2. We can also identify which rows and columns represent service sectors, as these tend to have a lower correlation and average influence values with non-service sectors (depicted in dark blue) and higher values among themselves.
The average adjusted mutual information obtained for the 1000 seed values tested is 0.90. This is a very high value which tells us that, on average, the partitions obtained for the directed and undirected versions of the network were very similar. This suggests that the community detection procedure is weakly dependent on the version of the network (directed vs. undirected) as well as the seed value used.

Discussion
In this paper we have introduced the biPCPG framework, a generalisation of the PCPG [30] algorithm to datasets with a multi-sample and multi-variable structure that allows a statistical significant and robust analysis, mainly by generating confidence bounds via an adapted bootstrapping procedure. We have then applied this new procedure to a recently introduced dataset that integrates the export of physical goods and services data. The proposed procedure allows the generation of a network of these economic sectors whose links represent the average influence in terms of temporal correlation. This can be seen as an an asymmetric formulation of relatedness [26,52]. The resulting network contains several hub nodes with high degree (namely Plastics, Pigments, Iron and steel articles, Preparations of cereals and milk and Aluminium) as well as distinct clusters of intuitively-related economic sectors (such as a food and plant cluster, a services cluster and manufacturing cluster). We find that, in this network, economic sectors display a relatively high assortativity according to their complexity rank and, to a lesser extent, their category.

Conclusions
In this work, we have introduced an asymmetric definition for relatedness by extending the PCPG methodology introduced in [30] for its use on bipartite datasets, which we call biPCPG. We apply this approach to a recently introduced dataset containing the exports of countries regarding both manufactured products and intangible services. We show that the biPCPG methodology is able to generate a statistically robust network of economic sectors which captures the underlying influence structure int erms of temporal correlations.
This work can be extended in a number of possible directions. First of all, the biPCPG framework can be applied to any temporal bipartite network, such as those of common use in economic complexity, such as the company-technology [9] or the country-scientific field network [29]. Moreover, the adapted bootstrapping procedure can be used to other network-generating techniques based on correlation-filtering to datasets with a multisample and multi-variable structure. These techniques include those based on threshold methods [53], the Minimum Spanning Tree [33] and the aforementioned PMFG [31], as well as more recent techniques based on a null-model approach [54]. This would be possible by replacing the last step in our procedure, the original PCPG algorithm, with the correlation-filtering technique of interest. Finally, it would also be particularly interesting to apply our procedure to datasets with the same structure but longer time series, such as financial datasets containing, for example, asset prices at the different exchanges where they are traded.  Acknowledgments: The authors acknowledge Michele Tumminello for providing the PCPG Mathematica code.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Fitness and Complexity of Economic sectors
From the matrices containing RCA c,p time series, described in Section 3.3 we can derive the M y matrix which has entries given by where c represents a country, p represents a product (or service), and y represents a given year. This matrix therefore summarises the countries having a comparative advantage at exporting the different products or services in a given year, or not. Two key quantities from the economic complexity literature are defined using this matrix, namely the fitness of countries and the complexity of products (or services) [17,55]. The intuition behind these quantities is that the higher the fitness of a country the higher its capability of exporting products of high complexity. It is therefore natural for the fitness to be proportional to the weighted sum of the products of which it is a competitive exporter. The definition of the complexity of a product is more subtle. In general terms, the complexity of a product should be inversely proportional to the number of countries exporting it. We should also note that more economically developed countries tend to have a highly diversified export basket, while less economically developed countries tend to have a much more limited diversification in their exports, and focused on low complexity products. Therefore, the upper bound of a product's complexity should be determined by the fitness of the countries' exporting it, with a strong bias towards lower fitness countries: if a product is exported by lower fitness countries, its complexity can not be high. The fitness F c of a country and the complexity Q p of a product (or service) are therefore defined using the following set of coupled iterative equations which are iterated until a fixed point is reached [56]. This fixed point has been shown to be stable and not dependent on the initial conditions, which are set toQ (0) p = 1∀p and F (0) c = 1∀c [17]. We use the complexity of products and services in our dataset to calculate an assortativity metric on the network G as described in Section 4.2.
It is worth noting that the dataset analysed and similar datasets explored in the economic complexity literature exhibit a nested structure [56]. This nested structure is manifested as a triangular structure in the M y matrices when countries (rows) and sectors (columns) are sorted by their fitness and complexity rank, respectively. This can be seen in Figure A1, which is the M y matrix for the year y = 2005.

Appendix B. Confidence and Prediction interval calculations
The 95% confidence interval around a linear fitμ y|x 0 done on n data points (x i , y i ) n = 1, . . . , n contains the mean response of new values µ y|x 0 at a given value x 0 with a 95% probability. This is given by whereμ y|x 0 = a + bx 0 is computed from the linear fit, T .975 n−2 is the 97.5th percentile of the Student's t-distribution with n − 2 degrees of freedom andσ is the standard deviation of the residuals in the linear fit given bŷ The 95% prediction interval around a linear fitŷ 0 is the interval within which a new observation, y 0 , at a given value, x 0 , is found, with 95% probability. This is given by whereŷ 0 = a + bx 0 is computed from the linear fit. See [57] for a more detailed description. Figure A2. Average correlation matrixK sorted by communities ν found by maximising modularity. Figure A3. Matrix showing average influence values between products d(p : p ) sorted by communities ν found by maximising modularity. Entries in white indicate that the average influence of a sector on itself is undefined.