Supporting Materials for the Manuscript:

A Markovian model of evolving world input-output network

Vahid Moosavi

Giulio Isacchini

ETH Zurich, Switzerland

(Under review at the PLOS ONE journal), 2017

For Direct Access to the Python notebooks and Data set, please go to the project repo

Manuscript can be found here

In [1]:
import datetime
import pandas as pd
# import pandas.io.data
import numpy as np
from matplotlib import pyplot as plt
import sys
import sompy as SOM# from pandas import Series, DataFrame
pd.__version__
%matplotlib inline
from pylab import matshow, savefig
from scipy.linalg import norm
import time

The proposed Markov model from input-output table

Here, we assume two economies each with one industry.
There is also one node corresponding to the Governments and Housholds within each economy

However, considering the changes in the values of Input-output models, we have a kind of time varying markov chain. Similarly, we assume
we have a time series of stationary markov chains, where for each year we fit one markov chain
In this experiment, we have focused on 41 economies (countries) with 35 industries within each economy interacting with other industries and economies
In [40]:
# Please, first run the Supporting Functions at the end of this notebook
In [9]:
economy_n = 41
industry_n = 35
header = pd.read_csv('./Data/WIOD/header.csv',header=None)
economy_names =header.values[1,range(0,1435,35)]
industry_names =header.values[0,range(0,35)]
Countries =pd.read_csv("./Data/WIOD/Economies_Names.csv")
# We calculated this before from the WIOD, for GDP calculations
VA = pd.read_csv('./Data/WIOD/VAs.csv',index_col=[0])
In [10]:
WIO = list()
for i in range(1995,2012):
    d = pd.read_csv('./Data/WIOD/wiot'+str(i)+'_row_apr12.csv',header=[0,1,2])
    WIO.append(d.values[:])
    print i, d.shape
1995 (1435, 1640)
1996 (1435, 1640)
1997 (1435, 1640)
1998 (1435, 1640)
1999 (1435, 1640)
2000 (1435, 1640)
2001 (1435, 1640)
2002 (1435, 1640)
2003 (1435, 1640)
2004 (1435, 1640)
2005 (1435, 1640)
2006 (1435, 1640)
2007 (1435, 1640)
2008 (1435, 1640)
2009 (1435, 1640)
2010 (1435, 1640)
2011 (1435, 1640)

Using Power Iteration to calculate the steady state probabilities, Pi

In [42]:
Pi =[]
TMs = []
Mixing_times =[]
singular_ids = []
fig3 = plt.figure(figsize=(15,7))
for i,WIOT in enumerate(WIO):
    economy_int_consumptions,industry_production_costs = aggregate_economy_int_consumptions(WIOT)
    economy_n = 41
    industry_n = 35
    states_n = economy_n*industry_n+economy_n
    TM = np.zeros((states_n,states_n))
    TM, singular_id = build_markov_chain_Dense(TM,WIOT,economy_n,industry_n,economy_int_consumptions,industry_production_costs)
    TMs.append(TM)
    singular_ids.extend(singular_id)
    t,mixing_ = simulate_markov(TM)
    Pi.append(t)
    Mixing_times.append(mixing_)

singular_ids = np.unique(singular_ids)
Making the TM in 0.592767953873 second
Iterating 290 times in 0.160306930542
Making the TM in 0.594426870346 second
Iterating 258 times in 0.145621061325
Making the TM in 0.583817958832 second
Iterating 250 times in 0.139415979385
Making the TM in 0.583045005798 second
Iterating 259 times in 0.143721103668
Making the TM in 0.581876993179 second
Iterating 261 times in 0.144848823547
Making the TM in 0.599870920181 second
Iterating 240 times in 0.136317014694
Making the TM in 0.592772960663 second
Iterating 228 times in 0.127722024918
Making the TM in 0.598742961884 second
Iterating 229 times in 0.127930879593
Making the TM in 0.586272954941 second
Iterating 221 times in 0.122807979584
Making the TM in 0.591900110245 second
Iterating 208 times in 0.115034818649
Making the TM in 0.590232133865 second
Iterating 190 times in 0.108383893967
Making the TM in 0.586740970612 second
Iterating 173 times in 0.0964498519897
Making the TM in 0.66259098053 second
Iterating 167 times in 0.0945069789886
Making the TM in 0.596711158752 second
Iterating 157 times in 0.0995650291443
Making the TM in 0.619015932083 second
Iterating 192 times in 0.108705043793
Making the TM in 0.608633041382 second
Iterating 173 times in 0.0964958667755
Making the TM in 0.583389043808 second
Iterating 157 times in 0.0878601074219
<matplotlib.figure.Figure at 0x13d027090>

Fig 2. The sequence of average mixing time of Markov chains as an aggregate index of globalization

Lower values indicate more globally connected network. The vertical bars show the error bars of minus/plus 3 standard deviations.

In [12]:
# We run the power iteration for several times to see if there are some variations in the mixing time
n =100
Mixing_times_n_times = np.zeros((n,len(TMs)))
for i in range(n):
    for j in range(len(TMs)):
        t,mixing_ = simulate_markov(TMs[j],verbose='off')
        Mixing_times_n_times[i,j]= mixing_
        
In [14]:
fig = plt.figure(figsize=(10,7))
yerr = 3*Mixing_times_n_times.std(axis=0)
plt.errorbar(range(1,18),Mixing_times_n_times.mean(axis=0), yerr=yerr,fmt='-',color='b',ecolor='b',linewidth=.51, elinewidth=3)
plt.xticks(range(1,18),range(1995,2012))
label = range(1995,2012)
font = {'size'   : 12}
plt.rc('font', **font)
plt.ylabel('Mixing Time of the Corresponding Markov Chain')
plt.xlabel('Year')
plt.xticks(range(1,18),label,rotation=90)

plt.grid(linewidth=.31,color='gray',linestyle='--')
path = './Images/forpaper/Mixing_Time.tiff'
# fig.savefig(path,dpi=300)
plt.show()

Fig 3. The sequence of Kemeny constants of Markov chains as an aggregate index of globalization

Lower values indicate more globally connected networks.

In [15]:
Kemenys = []
i = 0
for TM in TMs:
    K, pi = Kemeny_constant(TM)
    Kemenys.append(K)
    i+=1
    print i
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
In [16]:
fig = plt.figure(figsize=(10,7))
plt.plot(Kemenys,'.-b',linewidth=.51)
label = range(1995,2012)
plt.xticks(range(17),label,rotation=90)
plt.grid(linewidth=.31,color='gray',linestyle='--')
plt.xlabel('Year')
plt.ylabel('Kemeny Constant of the Corresponding Markov Chain')
font = {'size'   : 12}
plt.rc('font', **font)
path = './Images/forpaper/Kemeny.tiff'
# fig.savefig(path,dpi=300)

Fig 4. GDP shares of economies (red line) compared with their aggregated structural powers (blue line) over time

In [43]:
Data_all = np.asarray(Pi).T

fig = plt.figure(figsize=(18,18))
font = {'size'   : 8}
plt.rc('font', **font)
GDP_shares = np.zeros((economy_n,17))
Economies_shares = np.zeros((economy_n,17))
for i in range(economy_n):
    ind_indusrty =   range(i*industry_n,(i+1)*industry_n)
    GDP_shares[i]=VA.values[:,ind_indusrty].sum(axis=1)/VA.values[:].sum(axis=1)
    
    plt.subplot(7,6,i+1)
    plt.title(Countries.Country[i])
    ind_indusrty = range(i*industry_n,(i+1)*industry_n)+ [industry_n*economy_n+i]
    Economies_shares[i]=Data_all[ind_indusrty].T.sum(axis=1)
    
    plt.plot(range(17),Economies_shares[i],'-b',linewidth=1.2,label='Aggregated Structural Power of Economy')
    plt.plot(range(17),GDP_shares[i],'-r',linewidth=1.2,label='GDP Share of Economy')
    y1 = Economies_shares[i]
    y2 = GDP_shares[i]
    x = range(17)
    
    
    
    plt.fill_between(x, y1, y2, where=y2 >= y1, facecolor='red',alpha=.5, interpolate=True)
    plt.fill_between(x, y1, y2, where= y2 <= y1, facecolor='blue', alpha=.5,interpolate=True)
    path = './Images/forpaper/GDP_Share_Pi.tiff'
    label = range(1995,2012)
    plt.xticks(range(17),label,rotation=90)
    

    plt.tight_layout()
    
    plt.grid(linewidth=.31,color='gray',linestyle='--')

plt.legend(bbox_to_anchor=(2.29,1.022))
# fig.savefig(path,dpi=300)
Out[43]:
<matplotlib.legend.Legend at 0x13d018c50>

Fig 5. Predicted trends of structural potential of different economies

In [19]:
fig = plt.figure(figsize=(12,8))



min_x = np.min(Economies_shares[:,:]/GDP_shares[:,:])
min_y = np.min(GDP_shares[:,:])

max_x = np.max(Economies_shares[:,:]/GDP_shares[:,:])
max_y = np.max(GDP_shares[:,:])
eps = .02

for c in range(economy_n):

    
    x = Economies_shares[c,:]
    color = Economies_shares[c,:]/GDP_shares[c,:]
    x = range(17)
    y = GDP_shares[c,:]
    y = Economies_shares[c,:]
    y = Economies_shares[c,:]/GDP_shares[c,:]
    
    min_x = np.min(Economies_shares[c,:])
    min_y = np.min(GDP_shares[c,:])
    min_min = np.min([min_x,min_y])

    max_x = np.max(Economies_shares[c,:])
    max_y = np.max(GDP_shares[c,:])
    max_max = np.max([max_x,max_y]) 
    eps = .02
    

    ax = plt.gca()
    pred_period = 5
    col = 0
    if Countries.Country[c] in ["Brazil","Russia","China","India"]:
        ax.plot(x,y,'-',linewidth=1.2)
        
        
        all_preds = []
        for span in range(3,7):
            series = y.copy()
            
            y_pred  =update_by_momentum(series,span,pred_period)
            all_preds.append(y_pred)
            y_pred = list(y[-1:]) + y_pred
            x_pred = x[-1:]+ range(17,17+pred_period)
            ax.plot(x_pred,y_pred,'-r',linewidth=0.2)
        y_pred = list(y[-1:]) + list(np.median(all_preds,axis=0))
        x_pred = x[-1:]+ range(17,17+pred_period)
        ax.plot(x_pred,y_pred,'--r',linewidth=1.2)
        
        
        
        ax.scatter(x,y,vmax=1.7,vmin=0.5,c=color,s=20,marker='.',edgecolor='None', cmap='RdYlBu' ,alpha=1)
        ax.annotate(Countries.Country[c], size=10,xy = (x[16],y[-1]+.002), xytext = (0, 0),
                 textcoords = 'offset points', ha = 'right', va = 'bottom')
    else:
        ax.plot(x,y,'gray',linewidth=.3)
    
    

    ax.plot([0,16+pred_period],[1,1],'--k',linewidth=.2)
    
    plt.ylim([.25,2])
    plt.xlim([0,16+pred_period])
    
    
    plt.ylabel('Aggregated Structural Power of Economy Over GDP Share of Economy\n (Structural Potential of Economy)')
    plt.xlabel('Year')
    plt.xticks(range(17+pred_period),range(1995,2012+pred_period), rotation='vertical')

    plt.grid(linewidth=.31,color='gray',linestyle='--')
    plt.tight_layout()
    font = {'size'   : 10}
    plt.rc('font', **font)
    path = './Images/forpaper/GDP_Share_Vs_Pi_Share_trend.tiff'

# fig.savefig(path,dpi=300)