<jupyter_start><jupyter_text>Visualising time-seriesWe visualise the time-series to get a better idea of the characteristics of our interlinkages.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pickle
import pandas as pd
from sklearn.preprocessing import scale
from sklearn.metrics.pairwise import pairwise_kernels, pairwise_distances
from sklearn.decomposition import PCA, KernelPCA
import seaborn as sns
import dcor
from scipy.stats import pearsonr
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, Matern, Exponentiation, DotProduct, WhiteKernel
from sklearn.kernel_ridge import KernelRidge
from sklearn.linear_model import LinearRegression
import numpy.polynomial.polynomial as poly
import warnings 
warnings.filterwarnings('ignore')
indicators_values_i = pickle.load(open('utils/data/indicators_values_i_up_wb.pkl', 'rb'))
targets_values_i = pickle.load(open('utils/data/targets_values_i_up_arr_wb.pkl', 'rb'))
goals_values_i = pickle.load(open('utils/data/goals_values_i_up_arr_wb.pkl', 'rb'))
dict_all = pickle.load(open('utils/data/dict_all_wb.pkl', 'rb'))
dict_all_std = pickle.load(open('utils/data/dict_all_wb_std.pkl', 'rb'))
dict_all_i = pickle.load(open('utils/data/dict_all_i_wb.pkl', 'rb'))
averages = pickle.load(open('distance_cor/goals/grouping_averages.pkl', 'rb'))
dict_t = pickle.load(open('utils/data/temp.pkl', 'rb'))
# read amended csv file
c = pd.read_csv('utils/countries_wb.csv', dtype=str, delimiter=';', header=None)
countries = list(c[0])
continents = pd.read_csv(r'utils/continents.csv')
continents.replace({"Democratic People's Republic of Korea": "Korea, Dem. People's Rep.", 'Gambia': 'Gambia, The', 'United Kingdom of Great Britain and Northern Ireland': 'United Kingdom', 'Congo': 'Congo, Rep.', 'Democratic Republic of the Congo': 'Congo, Dem. Rep.', 'Czechia': 'Czech Republic', 'Iran (Islamic Republic of)': 'Iran, Islamic Rep.', "Côte d'Ivoire": "Cote d'Ivoire", 'Kyrgyzstan': 'Kyrgyz Republic', "Lao People's Democratic Republic": 'Lao PDR', 'Republic of Moldova': 'Moldova', 'Micronesia (Federated States of)': 'Micronesia, Fed. Sts.', 'Slovakia': 'Slovak Republic', 'Viet Nam': 'Vietnam', 'Egypt': 'Egypt, Arab Rep.', 'United Republic of Tanzania': 'Tanzania','United States of America': 'United States', 'Venezuela (Bolivarian Republic of)': 'Venezuela, RB', 'Yemen': 'Yemen, Rep.', 'Bahamas': 'Bahamas, The', 'Bolivia (Plurinational State of)': 'Bolivia'}, inplace=True)
groups = pd.read_csv(r'utils/groups.csv')
groups.replace({"Democratic People's Republic of Korea": "Korea, Dem. People's Rep.", 'Gambia': 'Gambia, The', 'United Kingdom of Great Britain and Northern Ireland': 'United Kingdom', 'Congo': 'Congo, Rep.', 'Democratic Republic of the Congo': 'Congo, Dem. Rep.', 'Czechia': 'Czech Republic', 'Iran (Islamic Republic of)': 'Iran, Islamic Rep.', "Côte d'Ivoire": "Cote d'Ivoire", 'Kyrgyzstan': 'Kyrgyz Republic', "Lao People's Democratic Republic": 'Lao PDR', 'Republic of Moldova': 'Moldova', 'Micronesia (Federated States of)': 'Micronesia, Fed. Sts.', 'Slovakia': 'Slovak Republic', 'Viet Nam': 'Vietnam', 'Egypt': 'Egypt, Arab Rep.', 'United Republic of Tanzania': 'Tanzania','United States of America': 'United States', 'Venezuela (Bolivarian Republic of)': 'Venezuela, RB', 'Yemen': 'Yemen, Rep.', 'Bahamas': 'Bahamas, The', 'Bolivia (Plurinational State of)': 'Bolivia'}, inplace=True)
info = pd.read_csv(r'utils/wb_info.csv', header=None)
period = ['2000', '2001', '2002', '2003', '2004', '2005', '2006', '2007', '2008', '2009', '2010', '2011', '2012', '2013', '2014', '2015', '2016', '2017', '2018', '2019']
info = info[info[0]!='DT.ODA.ODAT.CD']
# removes key in-place
countries.remove('Micronesia, Fed. Sts.')
continents['Oceania (excl. AUS + NZ)'] = continents['Oceania (excl. AUS + NZ)'].drop(index=4)
continents['Oceania (incl. AUS + NZ)'] = continents['Oceania (incl. AUS + NZ)'].drop(index=6)
continents['World'] = continents['World'].drop(index=170)
groups['Global South'].drop(index=1, inplace=True)
for country in countries:
    dict_t[country].set_index('YEAR', inplace=True)
    dict_t[country].index = dict_t[country].index.astype(str)
    dict_all_std[country] = dict_all_std[country].append(dict_t[country].AVG[-len(period):])
    dict_all_std[country].rename({'AVG': 'Temperature'}, inplace=True)
    dict_all[country] = dict_all[country].append(dict_t[country].AVG[-len(period):])
    dict_all[country].rename({'AVG': 'Temperature'}, inplace=True)
indicators = list(dict_all['France'].index)
targets = list(info[4].unique())

dict_targets = {}

for target in targets:
    t = info[0].where(info[4] == target)

    dict_targets[target] = [i for i in t if str(i) != 'nan']
#check 
dict_targets['7.1']
goals = list(info[3].unique())

dict_goals = {}

for goal in goals:
    g = info[4].where(info[3] == goal)

    dict_goals[goal] = [t for t in g if str(t) != 'nan']
    dict_goals[goal] = list(set(dict_goals[goal]))<jupyter_output><empty_output><jupyter_text>Defining combinations of goals:<jupyter_code>def combinations(iterable, r):
    # combinations('ABCD', 2) --> AB AC AD BC BD CD
    # combinations(range(4), 3) --> 012 013 023 123
    pool = tuple(iterable)
    n = len(pool)
    if r > n:
        return
    indices = list(range(r))
    yield list(pool[i] for i in indices)
    while True:
        for i in reversed(range(r)):
            if indices[i] != i + n - r:
                break
        else:
            return
        indices[i] += 1
        for j in range(i+1, r):
            indices[j] = indices[j-1] + 1
        yield list(pool[i] for i in indices)
g_combinations = list(combinations(goals, 2))
pairs_g1 = pd.DataFrame.from_records(data=g_combinations, columns=['pair_0', 'pair_1'])
# data preparation for continents (goals)
continents_prep_g = {}

for continent in continents:
    print(continent)
    
    continents_prep_g[continent] = np.empty(18, dtype=object)
    
    for g, goal in enumerate(goals):
        g_list = []
        for country in continents[continent].dropna():
            g_list.append(np.asarray(goals_values_i[country][g]))

        continents_prep_g[continent][g] = np.asarray(g_list)
# data preparation for continents (indicators)
continents_prep_i = {}

for continent in continents:
    print(continent)
    
    continents_prep_i[continent] = np.empty(len(dict_all['France']), dtype=object)
    
    for i, ind in enumerate(list(dict_all['France'].index)):
        i_list = []
        for country in continents[continent].dropna():
            i_list.append(np.asarray(dict_all[country].iloc[i].dropna()))

        continents_prep_i[continent][i] = np.asarray(i_list)
# data preparation for groups (goals)
groups_prep_g = {}

for group in groups:
    print(group)
    
    groups_prep_g[group] = np.empty(18, dtype=object)
    
    for g, goal in enumerate(goals):
        g_list = []
        for country in groups[group].dropna():
            g_list.append(np.asarray(goals_values_i[country][g]))

        groups_prep_g[group][g] = np.asarray(g_list)
# data preparation for groups (indicators)
groups_prep_i = {}

for group in groups:
    print(group)
    
    groups_prep_i[group] = np.empty(len(dict_all['France']), dtype=object)
    
    for i, ind in enumerate(list(dict_all['France'].index)):
        i_list = []
        for country in groups[group].dropna():
            i_list.append(np.asarray(dict_all[country].iloc[i].dropna()))

        groups_prep_i[group][i] = np.asarray(i_list)
# confidence intervals
def confidence(mean):
    interval = np.sqrt((mean * (100 - mean)) / 100)
    return interval

def width(Z):
    """
    Computes the median heuristic for the kernel bandwidth
    """
    dist_mat = pairwise_distances(Z, metric='euclidean')
    width_Z = np.median(dist_mat[dist_mat > 0])
    return width_Z<jupyter_output><empty_output><jupyter_text>Visualising dependence between nodes without the influence of the conditional setsTo start with, we visualise the dependence between two indicators of the variables $X$ and $Y$ with standardised and imputed data we have available for them. We take the latest available measurements per country.Our goal is to visualise the relationship between any two distance correlated nodes without the influence of the conditional set. We do so by regressing the conditional set $Z$ on both variables $X$ and $Y$ (with kernel ridge regression), and computing the residuals $r_x = x - f(z)$, and $r_y = y - f(z)$.Then, we take the residuals of the latest measurements of the two indicators of the variables $X$ and $Y$ and plot these again.<jupyter_code># conditions
HI_cond = pd.read_csv('distance_cor/goals/conditions_High Income.csv', index_col=[0])
UMI_cond = pd.read_csv('distance_cor/goals/conditions_Upper middle Income.csv', index_col=[0])
LMI_cond = pd.read_csv('distance_cor/goals/conditions_Lower middle Income.csv', index_col=[0])
LI_cond = pd.read_csv('distance_cor/goals/conditions_Low Income.csv', index_col=[0])
LAC_cond = pd.read_csv('distance_cor/goals/conditions_Latin America and the Caribbean.csv', index_col=[0])
OPEC_cond = pd.read_csv('distance_cor/goals/conditions_OPEC.csv', index_col=[0])
SE_cond = pd.read_csv('distance_cor/goals/conditions_Southern Europe.csv', index_col=[0])
GS_cond = pd.read_csv('distance_cor/goals/conditions_Global South.csv', index_col=[0])
GN_cond = pd.read_csv('distance_cor/goals/conditions_Global North.csv', index_col=[0])
# example: SDG 8 and 17 in High Income
GS_8_17 = GS_cond.loc[(GS_cond['pair_0']=='8') & (GS_cond['pair_1']=='17')]
GS_8_17<jupyter_output><empty_output><jupyter_text>Comparison of dependence before and after regressed on conditional subset  Example 1<jupyter_code># before
ind_0 = 'NY.GDP.PCAP.PP.KD' 
ind_1 = 'SE.PRM.UNER.ZS' 
group = 'Latin America and the Caribbean'
last_n = 0

x = []
y = []
for country in continents[group].dropna():
#for country in groups[group].dropna():
    if len(dict_all_i[country].loc[ind_0].dropna())>last_n and len(dict_all_i[country].loc[ind_1].dropna())>last_n:
        x.extend(dict_all_i[country].loc[ind_0].dropna().tail(last_n+1).to_list())
        y.extend(dict_all_i[country].loc[ind_1].dropna().tail(last_n+1).to_list())

if len(x)==0 or len(y)==0:
    pass
else:
    x = np.asarray(x).reshape(-1,1)
    y = np.asarray(y).reshape(-1,1)
    
    d_cor = dcor.distance_correlation(x, y)
    d_cor_test = dcor.independence.distance_covariance_test(x, y, num_resamples=10000)

    p_cor = pearsonr(x.reshape(-1), y.reshape(-1))

    # plot x against y; each country is one data point
    plt.figure(figsize=(18,14))
    plt.scatter(x, y, c='k')
    
    # fitting non-linear function to data
    kernel = Matern(nu=2.5)
    #kernel = RBF(length_scale=1)
    gpr_x = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=100, alpha=0.025)
    gpr_x.fit(x, y)
    x_ = np.linspace(min(x), max(x), 100)
    y_mean_gpr, y_std_gpr = gpr_x.predict(x_, return_std=True)
    
    plt.plot(x_, y_mean_gpr, '#adb0b3', lw=3, ls='-', label='Nonlinear fit')
    
    # fitting linear function to data
    linreg = LinearRegression()
    linreg.fit(x, y)
    y_mean_linreg = linreg.predict(x_)
    plt.plot(x_, y_mean_linreg, '#adb0b3', lw=3, ls='--', label='Linear fit')
    
    plt.xlabel(r'{}, standardised'.format(info[info[0]==ind_0][2].to_list()[0]), fontsize=24, labelpad=10)
    plt.ylabel(r'{}, standardised'.format(info[info[0]==ind_1][2].to_list()[0]), fontsize=24, labelpad=10)
    plt.xticks(size=18, alpha=0.65)
    plt.yticks(size=18, alpha=0.65)
    plt.title('{}'.format(group), y=1.015, fontdict={'fontsize': 44})
    
    textstr = '\n'.join((
        'Distance correlation = %.2f ($p$ = %.2f)' % (d_cor, d_cor_test.p_value),
        'Pearson correlation = %.2f ($p$ = %.2f)' % (p_cor[0], p_cor[1])))
    props = dict(boxstyle='round', facecolor='white', edgecolor='grey')
    plt.text(-0.85, -1.275, textstr, fontsize=26, bbox=props)
    
    plt.grid(alpha=0.35)
    plt.legend(loc=4, fontsize=22, edgecolor='grey')
    plt.tight_layout()
    plt.savefig('distance_cor/conditioning_before.png', format='png')
    plt.show()
# after
ind_0 = 'NY.GDP.PCAP.PP.KD'
pair_0 = int(info[info[0]==ind_0][3])
ind_1 = 'SE.PRM.UNER.ZS' 
pair_1 = int(info[info[0]==ind_1][3])
group = 'Latin America and the Caribbean'
last_n = 0

LAC_cond['pair_0'] = LAC_cond['pair_0'].astype(str)
if pair_0 > pair_1:
    group_cond = LAC_cond.loc[(LAC_cond['pair_0']=='{}'.format(pair_1)) & (LAC_cond['pair_1']=='{}'.format(pair_0))]
else:
    group_cond = LAC_cond.loc[(LAC_cond['pair_0']=='{}'.format(pair_0)) & (LAC_cond['pair_1']=='{}'.format(pair_1))]

if pair_0=='T':
    #x_res = groups_prep_g[group][17]
    x_res = continents_prep_g[group][17]
else:
    #x_res = groups_prep_g[group][int(pair_0)-1]
    x_res = continents_prep_g[group][int(pair_0)-1]
if pair_1=='T':
    #y_res = groups_prep_g[group][17]
    y_res = continents_prep_g[group][17]
else:
    #y_res = groups_prep_g[group][int(pair_1)-1]
    y_res = continents_prep_g[group][int(pair_1)-1]
    
# extracting conditional variables from column 'condition'
cond_list = []
for i in group_cond.condition.to_list()[0].split():
    newstr = ''.join((ch if ch in '0123456789.-eT' else ' ') for ch in i)
    cond_list.extend([i for i in newstr.split()])

condition = []
for c in cond_list:
    if c=='T':
        #condition.extend(groups_prep_g[group][17].T)
        condition.extend(continents_prep_g[group][17].T)
    else:
        #condition.extend(groups_prep_g[group][int(c)-1].T)
        condition.extend(continents_prep_g[group][int(c)-1].T)

z = np.asarray(condition).T

krr_x = KernelRidge(alpha=1, kernel='rbf', gamma=0.5/(width(z)**2))
krr_y = KernelRidge(alpha=1, kernel='rbf', gamma=0.5/(width(z)**2))
krr_x.fit(z, x_res)
krr_y.fit(z, y_res)

# project z on x and y
z_proj_x_res = krr_x.predict(z)
z_proj_y_res = krr_y.predict(z)

# residuals
r_x = x_res - z_proj_x_res
r_y = y_res - z_proj_y_res

x = []
y = []
for c, country in enumerate(continents[group].dropna()):
#for country in groups[group].dropna():
    k_x = np.where(x_res[c]==indicators_values_i[country].loc[ind_0].tail(last_n+1)[0])[0]
    k_y = np.where(y_res[c]==indicators_values_i[country].loc[ind_1].tail(last_n+1)[0])[0]
    x.extend(r_x[c][k_x])
    y.extend(r_y[c][k_y])

if len(x)==0 or len(y)==0:
    pass
else:
    x = np.asarray(x).reshape(-1,1)
    y = np.asarray(y).reshape(-1,1)
    
    p_cor = pearsonr(x.reshape(-1), y.reshape(-1))
    
    # dcor with entire series
    pair_0_array = continents_prep_g[group][pair_0-1]
    pair_1_array = continents_prep_g[group][pair_1-1]
    d_cor = dcor.partial_distance_correlation(pair_0_array, pair_1_array, z)
    d_cor_test = dcor.independence.partial_distance_covariance_test(pair_0_array, pair_1_array, z, num_resamples=10000)

    # plot x against y; each country is one data point
    plt.figure(figsize=(18,14))
    plt.scatter(x, y, c='k')
    
    # fitting non-linear function to data
    kernel = Matern(nu=2.5)
    #kernel = RBF(length_scale=1)
    gpr_x = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=100, alpha=0.025)
    gpr_x.fit(x, y)
    x_ = np.linspace(min(x), max(x), 100)
    y_mean_gpr, y_std_gpr = gpr_x.predict(x_, return_std=True)
    
    plt.plot(x_, y_mean_gpr, '#adb0b3', lw=3, ls='-', label='Nonlinear fit')
    
    # fitting linear function to data
    linreg = LinearRegression()
    linreg.fit(x, y)
    y_mean_linreg = linreg.predict(x_)
    plt.plot(x_, y_mean_linreg, '#adb0b3', lw=3, ls='--', label='Linear fit')
    
    plt.xlabel(r'Residuals of {}, standardised'.format(info[info[0]==ind_0][2].to_list()[0]), fontsize=24, labelpad=10)
    plt.ylabel(r'Residuals of {}, standardised'.format(info[info[0]==ind_1][2].to_list()[0]), fontsize=24, labelpad=10)
    plt.xticks(size=18, alpha=0.65)
    plt.yticks(size=18, alpha=0.65)
    plt.title('{}'.format(group), y=1.015, fontdict={'fontsize': 44})
    
    textstr = '\n'.join((
        'Distance correlation = %.2f ($p$ = %.2f)' % (d_cor, d_cor_test.p_value),
        'Pearson correlation = %.2f ($p$ = %.2f)' % (p_cor[0], p_cor[1])))
    props = dict(boxstyle='round', facecolor='white', edgecolor='grey')
    plt.text(-1.35, 0.31, textstr, fontsize=26, bbox=props)
    
    plt.grid(alpha=0.35)
    plt.legend(loc=4, fontsize=22, edgecolor='grey')
    plt.tight_layout()
    plt.savefig('distance_cor/conditioning_after.png', format='png')
    plt.show()<jupyter_output><empty_output><jupyter_text>Example 2<jupyter_code>np.where(x_res[c]==indicators_values_i[country].loc[ind_0].dropna().tail(last_n+1)[0])[0]
# before
ind_0 = 'SP.URB.TOTL.IN.ZS'
ind_1 = 'per_si_allsi.cov_pop_tot'
pair_0 = info[info[0]==ind_0][3].to_list()[0]
pair_1 = info[info[0]==ind_1][3].to_list()[0]
group = 'OPEC'
last_n = 0

x = []
y = []
#for country in continents[group].dropna():
for country in groups[group].dropna():
    if len(dict_all_i[country].loc[ind_0].dropna())>last_n and len(dict_all_i[country].loc[ind_1].dropna())>last_n:
        x.append(dict_all_i[country].loc[ind_0].dropna().to_list())
        y.append(dict_all_i[country].loc[ind_1].dropna().to_list())

x_arr = np.asarray(x)
y_arr = np.asarray(y)

d_cor = dcor.distance_correlation(x_arr, y_arr)
d_cor_test = dcor.independence.distance_covariance_test(x_arr, x_arr, num_resamples=10000)

# plot
plt.figure(figsize=(18,14))

dimensions = min(x_arr.shape[1], y_arr.shape[1])
for c, dim in enumerate(range(dimensions)):
    x_dim = x_arr[:, dim]
    y_dim = y_arr[:, dim]
    
    plt.scatter(x_dim, y_dim, color=sns.color_palette('Reds', dimensions+10)[c+10])

plt.xlabel(r'{}, standardised'.format(info[info[0]==ind_0][2].to_list()[0]), fontsize=24, labelpad=10)
plt.ylabel(r'{}, standardised'.format(info[info[0]==ind_1][2].to_list()[0]), fontsize=24, labelpad=10)
plt.xticks(size=18, alpha=0.65)
plt.yticks(size=18, alpha=0.65)
plt.title('{}'.format(group), y=1.015, fontdict={'fontsize': 44})

textstr = 'Distance correlation = %.2f ($p$ = %.2f)' % (d_cor, d_cor_test.p_value)
props = dict(boxstyle='round', facecolor='white', edgecolor='grey')
plt.text(-1.15, 3, textstr, fontsize=26, bbox=props)

plt.grid(alpha=0.335)
plt.tight_layout()
#plt.savefig('distance_cor/conditioning_before_2.png', format='png')
plt.show()
# after

OPEC_cond['pair_0'] = OPEC_cond['pair_0'].astype(str)
if pair_0 > pair_1:
    group_cond = OPEC_cond.loc[(OPEC_cond['pair_0']=='{}'.format(pair_1)) & (OPEC_cond['pair_1']=='{}'.format(pair_0))]
else:
    group_cond = OPEC_cond.loc[(OPEC_cond['pair_0']=='{}'.format(pair_0)) & (OPEC_cond['pair_1']=='{}'.format(pair_1))]

    
if pair_0=='T':
    x_res = groups_prep_g[group][17]
    #x_res = continents_prep_g[group][17]
else:
    x_res = groups_prep_g[group][int(pair_0)-1]
    #x_res = continents_prep_g[group][int(pair_0)-1]
if pair_1=='T':
    y_res = groups_prep_g[group][17]
    #y_res = continents_prep_g[group][17]
else:
    y_res = groups_prep_g[group][int(pair_1)-1]
    #y_res = continents_prep_g[group][int(pair_1)-1]

# extracting conditional variables from column 'condition'
cond_list = []
for i in group_cond.condition.to_list()[0].split():
    newstr = ''.join((ch if ch in '0123456789.-eT' else ' ') for ch in i)
    cond_list.extend([i for i in newstr.split()])

condition = []
for c in cond_list:
    if c=='T':
        condition.extend(groups_prep_g[group][17].T)
        #condition.extend(continents_prep_g[group][17].T)
    else:
        condition.extend(groups_prep_g[group][int(c)-1].T)
        #condition.extend(continents_prep_g[group][int(c)-1].T)

z = np.asarray(condition).T

krr_x = KernelRidge(alpha=1, kernel='rbf', gamma=0.5/(width(z)**2))
krr_y = KernelRidge(alpha=1, kernel='rbf', gamma=0.5/(width(z)**2))
krr_x.fit(z, x_res)
krr_y.fit(z, y_res)

# project z on x and y
z_proj_x_res = krr_x.predict(z)
z_proj_y_res = krr_y.predict(z)

# residuals
r_x = x_res - z_proj_x_res
r_y = y_res - z_proj_y_res

x = []
y = []
#for c, country in enumerate(continents[group].dropna()):
for c, country in enumerate(groups[group].dropna()):
    k_x = np.where(x_res[c]==indicators_values_i[country].loc[ind_0].tail(last_n+1)[0])[0]
    k_y = np.where(y_res[c]==indicators_values_i[country].loc[ind_1].tail(last_n+1)[0])[0]
    if len(k_x) > 1:
        k_x = [k_x[-1]]
    if len(k_y) > 1:
        k_y = [k_y[-1]]
    x.extend(r_x[c][k_x])
    y.extend(r_y[c][k_y])

if len(x)==0 or len(y)==0:
    pass
else:
    x = np.asarray(x).reshape(-1,1)
    y = np.asarray(y).reshape(-1,1)
    
    p_cor = pearsonr(x.reshape(-1), y.reshape(-1))
    
    # dcor with entire series
    #pair_0_array = continents_prep_g[group][pair_0-1]
    #pair_1_array = continents_prep_g[group][pair_1-1]
    pair_0_array = groups_prep_g[group][pair_0-1]
    pair_1_array = groups_prep_g[group][pair_1-1]
    d_cor = dcor.partial_distance_correlation(pair_0_array, pair_1_array, z)
    d_cor_test = dcor.independence.partial_distance_covariance_test(pair_0_array, pair_1_array, z, num_resamples=10000)

    # plot x against y; each country is one data point
    plt.figure(figsize=(18,14))
    plt.scatter(x, y, c='k')
    
    # fitting non-linear function to data
    kernel = Matern(nu=2.5)
    #kernel = RBF(length_scale=1)
    gpr_x = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=100, alpha=1e-2)
    gpr_x.fit(x, y)
    x_ = np.linspace(min(x), max(x), 100)
    y_mean_gpr, y_std_gpr = gpr_x.predict(x_, return_std=True)
    
    plt.plot(x_, y_mean_gpr, '#adb0b3', lw=3, ls='-', label='Nonlinear fit')
    plt.fill_between(x_.flatten(), y_mean_gpr.flatten()-1.96*y_std_gpr, y_mean_gpr.flatten()+1.96*y_std_gpr, interpolate=True, alpha=0.25, color='#adb0b3')
    
    # fitting linear function to data
    #linreg = LinearRegression()
    #linreg.fit(x, y)
    #y_mean_linreg = linreg.predict(x_)
    #plt.plot(x_, y_mean_linreg, '#adb0b3', lw=3, ls='--', label='Linear fit')
    
    plt.xlabel(r'Residuals of {}, standardised'.format(info[info[0]==ind_0][2].to_list()[0]), fontsize=24, labelpad=10)
    plt.ylabel('Residuals of\n{}, standardised'.format(info[info[0]==ind_1][2].to_list()[0]), fontsize=24, labelpad=10)
    plt.xticks(size=18, alpha=0.65)
    plt.yticks(size=18, alpha=0.65)
    plt.title('{}'.format(group), y=1.015, fontdict={'fontsize': 44})
    
    textstr = '\n'.join((
        'Distance correlation = %.2f ($p$ = %.2f)' % (d_cor, d_cor_test.p_value),
        'Pearson correlation = %.2f ($p$ = %.2f)' % (p_cor[0], p_cor[1])))
    textstr = 'Distance correlation = %.2f ($p$ = %.2f)' % (d_cor, d_cor_test.p_value)
    props = dict(boxstyle='round', facecolor='white', edgecolor='grey')
    plt.text(-0.415, -0.2035, textstr, fontsize=26, bbox=props)
    
    plt.grid(alpha=0.35)
    plt.legend(loc=4, fontsize=22, edgecolor='grey')
    plt.tight_layout()
    plt.savefig('distance_cor/conditioning_after_2.png', format='png')
    plt.show()<jupyter_output><empty_output><jupyter_text>Example 3<jupyter_code># before
ind_0 = 'SE.TER.ENRR'
ind_1 = 'ER.FSH.PROD.MT'
group = 'Southern Europe'
last_n = 0

x = []
y = []
for country in continents[group].dropna():
#for country in groups[group].dropna():
    if len(dict_all_i[country].loc[ind_0].dropna())>last_n and len(dict_all_i[country].loc[ind_1].dropna())>last_n:
        x.extend(dict_all_i[country].loc[ind_0].dropna().tail(last_n+1).to_list())
        y.extend(dict_all_i[country].loc[ind_1].dropna().tail(last_n+1).to_list())

if len(x)==0 or len(y)==0:
    pass
else:
    x = np.asarray(x).reshape(-1,1)
    y = np.asarray(y).reshape(-1,1)
    
    d_cor = dcor.distance_correlation(x, y)
    d_cor_test = dcor.independence.distance_covariance_test(x, y, num_resamples=10000)

    p_cor = pearsonr(x.reshape(-1), y.reshape(-1))

    # plot x against y; each country is one data point
    plt.figure(figsize=(18,14))
    plt.scatter(x, y, c='k')
    
    # fitting non-linear function to data
    kernel = Matern(nu=2.5)
    #kernel = RBF(length_scale=1)
    gpr_x = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=100, alpha=1)
    gpr_x.fit(x, y)
    x_ = np.linspace(min(x), max(x), 100)
    y_mean_gpr, y_std_gpr = gpr_x.predict(x_, return_std=True)
    
    plt.plot(x_, y_mean_gpr, '#adb0b3', lw=3, ls='-', label='Nonlinear fit')
    plt.fill_between(x_.flatten(), y_mean_gpr.flatten()-1.96*y_std_gpr, y_mean_gpr.flatten()+1.96*y_std_gpr, interpolate=True, alpha=0.25, color='#adb0b3')
        
    # fitting linear function to data
    linreg = LinearRegression()
    linreg.fit(x, y)
    y_mean_linreg = linreg.predict(x_)
    plt.plot(x_, y_mean_linreg, '#adb0b3', lw=3, ls='--', label='Linear fit')
    
    plt.xlabel(r'{}, standardised'.format(info[info[0]==ind_0][2].to_list()[0]), fontsize=24, labelpad=10)
    plt.ylabel(r'{}, standardised'.format(info[info[0]==ind_1][2].to_list()[0]), fontsize=24, labelpad=10)
    plt.xticks(size=18, alpha=0.65)
    plt.yticks(size=18, alpha=0.65)
    plt.title('{}'.format(group), y=1.015, fontdict={'fontsize': 44})
    
    textstr = '\n'.join((
        'Distance correlation = %.2f ($p$ = %.2f)' % (d_cor, d_cor_test.p_value),
        'Pearson correlation = %.2f ($p$ = %.2f)' % (p_cor[0], p_cor[1])))
    props = dict(boxstyle='round', facecolor='white', edgecolor='grey')
    plt.text(4.1, -1.5, textstr, fontsize=26, bbox=props)
    
    plt.grid(alpha=0.335)
    plt.legend(loc=4, fontsize=22, edgecolor='grey')
    plt.tight_layout()
    plt.savefig('distance_cor/conditioning_before_3.png', format='png')
    plt.show()
# after
pair_0 = int(info[info[0]==ind_0][3])
pair_1 = int(info[info[0]==ind_1][3])
last_n = 0

SE_cond['pair_0'] = SE_cond['pair_0'].astype(str)
group_cond = SE_cond.loc[(SE_cond['pair_0']=='{}'.format(pair_0)) & (SE_cond['pair_1']=='{}'.format(pair_1))]
   
if pair_0=='T':
    #x_res = groups_prep_g[group][17]
    x_res = continents_prep_g[group][17]
else:
    #x_res = groups_prep_g[group][int(pair_0)-1]
    x_res = continents_prep_g[group][int(pair_0)-1]
if pair_1=='T':
    pair_1 = 18
    #y_res = groups_prep_g[group][17]
    y_res = continents_prep_g[group][17]
else:
    #y_res = groups_prep_g[group][int(pair_1)-1]
    y_res = continents_prep_g[group][int(pair_1)-1]

# extracting conditional variables from column 'condition'
cond_list = []
for i in group_cond.condition.to_list()[0].split():
    newstr = ''.join((ch if ch in '0123456789.-eT' else ' ') for ch in i)
    cond_list.extend([i for i in newstr.split()])

condition = []
for c in cond_list:
    if c=='T':
        #condition.extend(groups_prep_g[group][17].T)
        condition.extend(continents_prep_g[group][17].T)
    else:
        #condition.extend(groups_prep_g[group][int(c)-1].T)
        condition.extend(continents_prep_g[group][int(c)-1].T)

z = np.asarray(condition).T
krr_x = KernelRidge(alpha=1, kernel='rbf', gamma=0.5/(width(z)**2))
krr_y = KernelRidge(alpha=1, kernel='rbf', gamma=0.5/(width(z)**2))
krr_x.fit(z, x_res)
krr_y.fit(z, y_res)

# project z on x and y
z_proj_x_res = krr_x.predict(z)
z_proj_y_res = krr_y.predict(z)

# residuals
r_x = x_res - z_proj_x_res
r_y = y_res - z_proj_y_res

x = []
y = []
for c, country in enumerate(continents[group].dropna()):
#for c, country in enumerate(groups[group].dropna()):
    k_x = np.where(x_res[c]==indicators_values_i[country].loc[ind_0].dropna().tail(last_n+1)[0])[0]
    k_y = np.where(y_res[c]==indicators_values_i[country].loc[ind_1].dropna().tail(last_n+1)[0])[0]
    if len(k_x) > 1:
        k_x = [k_x[-1]]
    if len(k_y) > 1:
        k_y = [k_y[-1]]
    x.extend(r_x[c][k_x])
    y.extend(r_y[c][k_y])

if len(x)==0 or len(y)==0:
    pass
else:
    x = np.asarray(x).reshape(-1,1)
    y = np.asarray(y).reshape(-1,1)
    
    p_cor = pearsonr(x.reshape(-1), y.reshape(-1))
    
    # dcor with entire series
    pair_0_array = continents_prep_g[group][pair_0-1]
    pair_1_array = continents_prep_g[group][pair_1-1]
    #pair_0_array = groups_prep_g[group][pair_0-1]
    #pair_1_array = groups_prep_g[group][pair_1-1]
    d_cor = dcor.partial_distance_correlation(pair_0_array, pair_1_array, z)
    d_cor_test = dcor.independence.partial_distance_covariance_test(pair_0_array, pair_1_array, z, num_resamples=10000)

    # plot x against y; each country is one data point
    plt.figure(figsize=(18,14))
    plt.scatter(x, y, c='k')
    
    # fitting non-linear function to data
    kernel = Matern(nu=2.5)
    #kernel = RBF(length_scale=1)
    gpr_x = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=100, alpha=1)
    gpr_x.fit(x, y)
    x_ = np.linspace(min(x), max(x), 100)
    y_mean_gpr, y_std_gpr = gpr_x.predict(x_, return_std=True)
    
    plt.plot(x_, y_mean_gpr, '#adb0b3', lw=3, ls='-', label='Nonlinear fit')
    plt.fill_between(x_.flatten(), y_mean_gpr.flatten()-1.96*y_std_gpr, y_mean_gpr.flatten()+1.96*y_std_gpr, interpolate=True, alpha=0.25, color='#adb0b3')
    
    # fitting linear function to data
    linreg = LinearRegression()
    linreg.fit(x, y)
    y_mean_linreg = linreg.predict(x_)
    plt.plot(x_, y_mean_linreg, '#adb0b3', lw=3, ls='--', label='Linear fit')
    
    plt.xlabel(r'Residuals of {}, standardised'.format(info[info[0]==ind_0][2].to_list()[0]), fontsize=24, labelpad=10)
    plt.ylabel('Residuals of {}, standardised'.format(info[info[0]==ind_1][2].to_list()[0]), fontsize=24, labelpad=10)
    plt.xticks(size=18, alpha=0.65)
    plt.yticks(size=18, alpha=0.65)
    plt.ylim((-1.75, 1.75))
    plt.title('{}'.format(group), y=1.015, fontdict={'fontsize': 44})
    
    textstr = '\n'.join((
        'Distance correlation = %.2f ($p$ = %.2f)' % (d_cor, d_cor_test.p_value),
        'Pearson correlation = %.2f ($p$ = %.2f)' % (p_cor[0], p_cor[1])))
    props = dict(boxstyle='round', facecolor='white', edgecolor='grey')
    plt.text(-0.46, -1.65, textstr, fontsize=26, bbox=props)
    
    plt.grid(alpha=0.35)
    plt.legend(loc=4, fontsize=22, edgecolor='grey')
    plt.tight_layout()
    plt.savefig('distance_cor/conditioning_after_3.png', format='png')
    plt.show()<jupyter_output><empty_output><jupyter_text>Comparing Global South and Global NorthWe plot the interlinkages particularly different in both groupings once the conditional sets are regressed out.<jupyter_code>info[info[3]=='10']
ind_1 = 'EN.ATM.CO2E.KD.GD'
ind_0 = 'EG.ELC.RNEW.ZS'
last_n = 0

pair_0 = int(info[info[0]==ind_0][3])
pair_1 = int(info[info[0]==ind_1][3])

GS_cond['pair_0'] = GS_cond['pair_0'].astype(str)
group_GS_cond = GS_cond.loc[(GS_cond['pair_0']=='{}'.format(pair_0)) & (GS_cond['pair_1']=='{}'.format(pair_1))]

GN_cond['pair_0'] = GN_cond['pair_0'].astype(str)
group_GN_cond = GN_cond.loc[(GN_cond['pair_0']=='{}'.format(pair_0)) & (GN_cond['pair_1']=='{}'.format(pair_1))]
   
if pair_0=='T':
    x_res_GS = groups_prep_g['Global South'][17]
    x_res_GN = groups_prep_g['Global North'][17]
    #x_res = continents_prep_g[group][17]
else:
    x_res_GS = groups_prep_g['Global South'][int(pair_0)-1]
    x_res_GN = groups_prep_g['Global North'][int(pair_0)-1]
    #x_res = continents_prep_g[group][int(pair_0)-1]
if pair_1=='T':
    pair_1 = 18
    y_res_GS = groups_prep_g['Global South'][17]
    y_res_GN = groups_prep_g['Global North'][17]
    #y_res = continents_prep_g[group][17]
else:
    y_res_GS = groups_prep_g['Global South'][int(pair_1)-1]
    y_res_GN = groups_prep_g['Global North'][int(pair_1)-1]
    #y_res = continents_prep_g[group][int(pair_1)-1]

# extracting conditional variables from column 'condition'
cond_list_GS = []
for i in group_GS_cond.condition.to_list()[0].split():
    newstr_GS = ''.join((ch if ch in '0123456789.-eT' else ' ') for ch in i)
    cond_list_GS.extend([i for i in newstr_GS.split()])

condition_GS = []
for c in cond_list_GS:
    if c=='T':
        condition_GS.extend(groups_prep_g['Global South'][17].T)
        #condition.extend(continents_prep_g[group][17].T)
    else:
        condition_GS.extend(groups_prep_g['Global South'][int(c)-1].T)
        #condition.extend(continents_prep_g[group][int(c)-1].T)

cond_list_GN = []
for i in group_GN_cond.condition.to_list()[0].split():
    newstr_GN = ''.join((ch if ch in '0123456789.-eT' else ' ') for ch in i)
    cond_list_GN.extend([i for i in newstr_GN.split()])

condition_GN = []
for c in cond_list_GN:
    if c=='T':
        condition_GN.extend(groups_prep_g['Global North'][17].T)
        #condition.extend(continents_prep_g[group][17].T)
    else:
        condition_GN.extend(groups_prep_g['Global North'][int(c)-1].T)
        #condition.extend(continents_prep_g[group][int(c)-1].T)

z_GS = np.asarray(condition_GS).T
krr_x_GS = KernelRidge(alpha=10, kernel='rbf', gamma=0.5/(width(z_GS)**2))
krr_y_GS = KernelRidge(alpha=10, kernel='rbf', gamma=0.5/(width(z_GS)**2))
krr_x_GS.fit(z_GS, x_res_GS)
krr_y_GS.fit(z_GS, y_res_GS)

z_GN = np.asarray(condition_GN).T
krr_x_GN = KernelRidge(alpha=10, kernel='rbf', gamma=0.5/(width(z_GN)**2))
krr_y_GN = KernelRidge(alpha=10, kernel='rbf', gamma=0.5/(width(z_GN)**2))
krr_x_GN.fit(z_GN, x_res_GN)
krr_y_GN.fit(z_GN, y_res_GN)

# project z on x and y
z_proj_x_res_GS = krr_x_GS.predict(z_GS)
z_proj_y_res_GS = krr_y_GS.predict(z_GS)

z_proj_x_res_GN = krr_x_GN.predict(z_GN)
z_proj_y_res_GN = krr_y_GN.predict(z_GN)

# residuals
r_x_GS = x_res_GS - z_proj_x_res_GS
r_y_GS = y_res_GS - z_proj_y_res_GS

r_x_GN = x_res_GN - z_proj_x_res_GN
r_y_GN = y_res_GN - z_proj_y_res_GN
 
# dcor with entire series
#pair_0_array = continents_prep_g[group][pair_0-1]
#pair_1_array = continents_prep_g[group][pair_1-1]
pair_0_array_GS = groups_prep_g['Global South'][pair_0-1]
pair_1_array_GS = groups_prep_g['Global South'][pair_1-1]
d_cor_GS = dcor.partial_distance_correlation(pair_0_array_GS, pair_1_array_GS, z_GS)
d_cor_test_GS = dcor.independence.partial_distance_covariance_test(pair_0_array_GS, pair_1_array_GS, z_GS, num_resamples=10000)

pair_0_array_GN = groups_prep_g['Global North'][pair_0-1]
pair_1_array_GN = groups_prep_g['Global North'][pair_1-1]
d_cor_GN = dcor.partial_distance_correlation(pair_0_array_GN, pair_1_array_GN, z_GN)
d_cor_test_GN = dcor.independence.partial_distance_covariance_test(pair_0_array_GN, pair_1_array_GN, z_GN, num_resamples=10000)

# plot x against y; each country is one data point
plt.figure(figsize=(18,14))

# finding indicator in pair_0_array again
len_pair_0 = 0
for ind in info[info[3]==str(pair_0)][0].tolist():
    l = len(indicators_values_i['Brazil'].loc[ind].dropna())
    len_pair_0 += l
    if ind==ind_0:
        plot_x_GS = r_x_GS[:, len_pair_0-l:len_pair_0]
        
len_pair_0 = 0
for ind in info[info[3]==str(pair_0)][0].tolist():
    l = len(indicators_values_i['France'].loc[ind].dropna())
    len_pair_0 += l
    if ind==ind_0:
        plot_x_GN = r_x_GN[:, len_pair_0-l:len_pair_0]

# finding indicator in pair_1_array again
len_pair_1 = 0
for ind in info[info[3]==str(pair_1)][0].tolist():
    l = len(indicators_values_i['Brazil'].loc[ind].dropna())
    len_pair_1 += l
    if ind==ind_1:
        plot_y_GS = r_y_GS[:, len_pair_1-l:len_pair_1]

len_pair_1 = 0
for ind in info[info[3]==str(pair_1)][0].tolist():
    l = len(indicators_values_i['France'].loc[ind].dropna())
    len_pair_1 += l
    if ind==ind_1:
        plot_y_GN = r_y_GN[:, len_pair_1-l:len_pair_1]
        
plt.scatter(plot_x_GS[:,-1], plot_y_GS[:,-1], c='forestgreen', s=60, label='Global South, distance correlation = %.2f ($p$ = %.2f)' % (d_cor_GS, d_cor_test_GS.p_value))
plt.scatter(plot_x_GN[:,-1], plot_y_GN[:,-1], c='orange', s=60, label='Global North, distance correlation = %.2f ($p$ = %.2f)' % (d_cor_GN, d_cor_test_GN.p_value))

# fitting non-linear function to data
#kernel = Matern(nu=2.5)
#kernel = RBF(length_scale=1)
#gpr_x = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=100, alpha=1)
#gpr_x.fit(x, y)
#x_ = np.linspace(min(x), max(x), 100)
#y_mean_gpr, y_std_gpr = gpr_x.predict(x_, return_std=True)

#plt.plot(x_, y_mean_gpr, '#adb0b3', lw=3, ls='-', label='Nonlinear fit')
#plt.fill_between(x_.flatten(), y_mean_gpr.flatten()-1.96*y_std_gpr, y_mean_gpr.flatten()+1.96*y_std_gpr, interpolate=True, alpha=0.25, color='#adb0b3')

#plt.xlabel(r'Residuals of {}, standardised'.format(info[info[0]==ind_0][2].to_GSst()[0]), fontsize=24, labelpad=10)
#plt.ylabel('Residuals of {}, standardised'.format(info[info[0]==ind_1][2].to_GSst()[0]), fontsize=24, labelpad=10)
plt.xticks(size=18, alpha=0.65)
plt.yticks(size=18, alpha=0.65)
plt.ylim((-1.75, 1.75))
plt.title('SDG {} vs SDG {}'.format(pair_0, pair_1), y=1.015, fontdict={'fontsize': 44})

plt.grid(alpha=0.35)
plt.legend(loc=4, fontsize=22, edgecolor='grey')
plt.tight_layout()
#plt.savefig('distance_cor/conditioning_GSGN_1.png', format='png')
plt.show()<jupyter_output><empty_output><jupyter_text>Plotting the first principle components of the residuals of $X$ and $Y$<jupyter_code>def residuals_kpca(pair_0, pair_1, group):
    
    if group=='High Income':
        HI_cond['pair_0'] = HI_cond['pair_0'].astype(str)
        group_cond = HI_cond.loc[(HI_cond['pair_0']=='{}'.format(pair_0)) & (HI_cond['pair_1']=='{}'.format(pair_1))]
    elif group=='Upper middle Income':
        UMI_cond['pair_0'] = UMI_cond['pair_0'].astype(str)
        group_cond = UMI_cond.loc[(UMI_cond['pair_0']=='{}'.format(pair_0)) & (UMI_cond['pair_1']=='{}'.format(pair_1))]
    elif group=='Lower middle Income':
        LMI_cond['pair_0'] = LMI_cond['pair_0'].astype(str)
        group_cond = LMI_cond.loc[(LMI_cond['pair_0']=='{}'.format(pair_0)) & (LMI_cond['pair_1']=='{}'.format(pair_1))]
    elif group=='Low Income':
        LI_cond['pair_0'] = LI_cond['pair_0'].astype(str)
        group_cond = LI_cond.loc[(LI_cond['pair_0']=='{}'.format(pair_0)) & (LI_cond['pair_1']=='{}'.format(pair_1))]
    elif group=='Latin America and the Caribbean':
        LAC_cond['pair_0'] = LAC_cond['pair_0'].astype(str)
        group_cond = LAC_cond.loc[(LAC_cond['pair_0']=='{}'.format(pair_0)) & (LAC_cond['pair_1']=='{}'.format(pair_1))]
    else:
        raise ValueError('Not possible for this group yet.')
    
    if pair_0=='T':
        #x = groups_prep_g[group][17]
        x = continents_prep_g[group][17]
    else:
        #x = groups_prep_g[group][int(pair_0)-1]
        x = continents_prep_g[group][int(pair_0)-1]
        
    if pair_1=='T':
        #y = groups_prep_g[group][17]
        y = continents_prep_g[group][17]
    else:
        #y = groups_prep_g[group][int(pair_1)-1]
        y = continents_prep_g[group][int(pair_0)-1]
    
    # extracting conditional variables from column 'condition'
    cond_list = []
    for i in group_cond.condition.to_list()[0].split():
        newstr = ''.join((ch if ch in '0123456789.-eT' else ' ') for ch in i)
        cond_list.extend([i for i in newstr.split()])

    condition = []
    for c in cond_list:
        if c=='T':
            #condition.extend(groups_prep_g[group][17].T)
            condition.extend(continents_prep_g[group][17].T)
        else:
            #condition.extend(groups_prep_g[group][int(c)-1].T)
            condition.extend(continents_prep_g[group][int(c)-1].T)

    z = np.asarray(condition).T

    krr_x = KernelRidge(alpha=1, kernel='rbf', gamma=0.5/(width(z)**2))
    krr_y = KernelRidge(alpha=1, kernel='rbf', gamma=0.5/(width(z)**2))
    krr_x.fit(z, x)
    krr_y.fit(z, y)
    z_proj_x = krr_x.predict(z)
    z_proj_y = krr_y.predict(z)

    # residuals
    r_x = x - z_proj_x
    r_y = y - z_proj_y 
    
    # selecting the principle component of r_y and r_x
    kpca_x = KernelPCA(n_components=1, kernel='rbf')
    kpca_y = KernelPCA(n_components=1, kernel='rbf')
    r_x_kpca = kpca_x.fit_transform(r_x)
    r_y_kpca = kpca_y.fit_transform(r_y)

    # plot r_x_kpca against r_y_kpca; each country is one data point
    plt.figure(figsize=(18,14))
    plt.scatter(r_x_kpca, r_y_kpca, c='k')

    kernel = Matern(nu=2.5)
    gpr_x = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=50, alpha=1e-2).fit(r_x_kpca, r_y_kpca)
    x_ = np.linspace(min(r_x_kpca), max(r_x_kpca), 100)
    y_mean, y_std = gpr_x.predict(x_, return_std=True)
    
    plt.plot(x_, y_mean, '#adb0b3', lw=3)
    plt.fill_between(x_.flatten(), y_mean.flatten()-1.96*y_std, y_mean.flatten()+1.96*y_std, interpolate=True, alpha=0.25, color='#adb0b3')
    x_label_x = 'SDG {}'.format(pair_0)
    if pair_1=='T':
        y_label_x = 'Mean Temperature'
    else:
        y_label_x = 'SDG {}'.format(pair_1)
    plt.xlabel(r'$1^{st}$ Principle Component of %s' % x_label_x, fontsize=32, labelpad=10)
    plt.ylabel(r'$1^{st}$ Principle Component of %s' % y_label_x, fontsize=32, labelpad=10)
    plt.xticks([])
    plt.yticks([])
    plt.title('{}'.format(group), y=1.015, fontdict={'fontsize': 44})
    plt.show()
    
    return r_x_kpca, r_y_kpca
# example
r_x_kpca, r_y_kpca = residuals_kpca(4, 17, 'Latin America and the Caribbean')<jupyter_output><empty_output><jupyter_text>Figure 1<jupyter_code>def corr(ind_0, ind_1, group, last_n=0):
    x = []
    y = []
    for country in continents[group].dropna():
    #for country in groups[group].dropna():
        if len(dict_all[country].loc[ind_0].dropna())>last_n and len(dict_all[country].loc[ind_1].dropna())>last_n:
            x.extend(dict_all[country].loc[ind_0].dropna().tail(last_n+1).to_list())
            y.extend(dict_all[country].loc[ind_1].dropna().tail(last_n+1).to_list())

    if len(x)<2 or len(y)<2:
        return None
    else:
        x = np.asarray(x).reshape(-1,1)
        y = np.asarray(y).reshape(-1,1)

        d_cor = dcor.distance_correlation(x, y)
        d_cor_test = dcor.independence.distance_covariance_test(x, y, num_resamples=5000)

        p_cor = pearsonr(x.reshape(-1), y.reshape(-1))
        
        return d_cor, d_cor_test.p_value, p_cor[0], p_cor[1]
examples = {}

for group in ['Latin America and the Caribbean']:  #groups.columns.to_list():
    print('----------------------------------------')
    print(group)
    examples[group] = {}
    
    for ind_0 in indicators:
        for ind_1 in indicators:
            if ind_0!=ind_1:
                c = corr(ind_0, ind_1, group)
                if c==None:
                    pass
                elif 0.4<c[0]<1 and c[1]<0.05 and -0.15<c[2]<0.15 and c[3]>0.4:
                    examples[group][(ind_0, ind_1)] = c
                    print(group, ind_0, ind_1, c)
                else:
                    pass
# save
exam = open('Figure_1_example_{}.pkl'.format(group), 'wb')
pickle.dump(examples[group], exam)
exam.close()
# load
exam = pickle.load(open('Figure_1_example_Latin America and the Caribbean.pkl', 'rb'))
ind_0 = 'NY.GDP.PCAP.PP.KD' 
ind_1 = 'SE.PRM.UNER.ZS' 
group = 'Latin America and the Caribbean'
last_n = 0

print(info[info[0]==ind_0][3])
print(info[info[0]==ind_1][3])

x = []
y = []
for country in continents[group].dropna():
#for country in groups[group].dropna():
    if len(dict_all[country].loc[ind_0].dropna())>last_n and len(dict_all[country].loc[ind_1].dropna())>last_n:
        x.extend(dict_all[country].loc[ind_0].dropna().tail(last_n+1).to_list())
        y.extend(dict_all[country].loc[ind_1].dropna().tail(last_n+1).to_list())

if len(x)==0 or len(y)==0:
    pass
else:
    x = np.asarray(x).reshape(-1,1)
    y = np.asarray(y).reshape(-1,1)
    
    d_cor = dcor.distance_correlation(x, y)
    d_cor_test = dcor.independence.distance_covariance_test(x, y, num_resamples=10000)

    p_cor = pearsonr(x.reshape(-1), y.reshape(-1))

    # plot x against y; each country is one data point
    plt.figure(figsize=(18,14))
    plt.scatter(x, y, c='k')
    
    # fitting non-linear function to data
    kernel = Matern(nu=2.5)
    #kernel = RBF(length_scale=1)
    gpr_x = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=100, alpha=0.435)
    gpr_x.fit(x, y)
    x_ = np.linspace(min(x), max(x), 100)
    y_mean_gpr, y_std_gpr = gpr_x.predict(x_, return_std=True)
    
    plt.plot(x_, y_mean_gpr, '#adb0b3', lw=3, ls='-', label='Nonlinear fit')
    #plt.fill_between(x_.flatten(), y_mean_gpr.flatten()-19.6*y_std_gpr.flatten(), y_mean_gpr.flatten()+19.6*y_std_gpr.flatten(), interpolate=True, alpha=0.25, color='#adb0b3')
    
    # fitting linear function to data
    linreg = LinearRegression()
    linreg.fit(x, y)
    y_mean_linreg = linreg.predict(x_)
    plt.plot(x_, y_mean_linreg, '#adb0b3', lw=3, ls='--', label='Linear fit')
    
    plt.xlabel(r'{}'.format(info[info[0]==ind_0][2].to_list()[0]), fontsize=24, labelpad=10)
    plt.ylabel(r'{}'.format(info[info[0]==ind_1][2].to_list()[0]), fontsize=24, labelpad=10)
    plt.xticks(size=18, alpha=0.65)
    plt.yticks(size=18, alpha=0.65)
    plt.title('{}'.format(group), y=1.015, fontdict={'fontsize': 44})
    
    textstr = '\n'.join((
        'Distance correlation = %.2f ($p$ = %.2f)' % (d_cor, d_cor_test.p_value),
        'Pearson correlation = %.2f ($p$ = %.2f)' % (p_cor[0], p_cor[1])))
    props = dict(boxstyle='round', facecolor='white', edgecolor='grey')
    plt.text(20700, 40, textstr, fontsize=26, bbox=props)
    
    #x_arrow, y_arrow = 42, 53
    #plt.arrow(x_arrow, y_arrow, 1, 6, head_width=0.75, head_length=0.75, fc='#1b75bc', ec='#1b75bc', width=0.2)
    #plt.text(x_arrow-1, y_arrow-3, 'T', fontsize=28, c='#1b75bc', fontfamily='serif', fontweight='roman')
    
    x_arrow_c1, y_arrow_c1 = 3400, 41
    plt.arrow(x_arrow_c1, y_arrow_c1, -800, +0.5, head_width=0.5, head_length=500, fc='orange', ec='orange', width=0.15)
    plt.text(x_arrow_c1+100, y_arrow_c1-0.5, 'Haiti', fontsize=28, c='orange', fontfamily='serif', fontweight='roman')
    
    x_arrow_c2, y_arrow_c2 = 34000, 24.5
    plt.arrow(x_arrow_c2, y_arrow_c2, 800, -0.5, head_width=0.5, head_length=500, fc='orange', ec='orange', width=0.15)
    plt.text(x_arrow_c2-6250, y_arrow_c2-0.25, 'The Bahamas', fontsize=28, c='orange', fontfamily='serif', fontweight='roman')
    
    x_arrow_c3, y_arrow_c3 = 33000, 20
    plt.arrow(x_arrow_c3, y_arrow_c3, 800, 0.5, head_width=0.5, head_length=500, fc='orange', ec='orange', width=0.15)
    plt.text(x_arrow_c3-5500, y_arrow_c3-0.5, '<NAME>', fontsize=28, c='orange', fontfamily='serif', fontweight='roman')
    
    x_arrow_c4, y_arrow_c4 = 20000, 3.25
    plt.arrow(x_arrow_c4, y_arrow_c4, 1000, -2, head_width=0.65, head_length=350, fc='orange', ec='orange', width=0.15)
    plt.text(x_arrow_c4-2000, y_arrow_c4+0.5, 'Argentina', fontsize=28, c='orange', fontfamily='serif', fontweight='roman')
    
    x_arrow_c5, y_arrow_c5 = 30000, 12.25
    plt.arrow(x_arrow_c5, y_arrow_c5, 800, 0.5, head_width=0.5, head_length=500, fc='orange', ec='orange', width=0.15)
    plt.text(x_arrow_c5-3600, y_arrow_c5-0.5, 'Panama', fontsize=28, c='orange', fontfamily='serif', fontweight='roman')
    
    plt.grid(alpha=0.35)
    plt.legend(loc=4, fontsize=22, edgecolor='grey')
    plt.tight_layout()
    plt.savefig('distance_cor/syn_tra_dependence_countries3.png', format='png')
    plt.show()
# countries at threshold
for country in continents[group].dropna():
    if dict_all[country].loc[ind_0].isnull().all()==False and dict_all[country].loc[ind_1].isnull().all()==False:
        print(country, dict_all[country].loc[ind_0][dict_all[country].loc[ind_0].last_valid_index()], dict_all[country].loc[ind_1][dict_all[country].loc[ind_1].last_valid_index()])<jupyter_output><empty_output><jupyter_text>Table with most important variables<jupyter_code>centr_continents = pickle.load(open('distance_cor/goals/centrality_continents.pkl', 'rb'))
centr_groups = pickle.load(open('distance_cor/goals/centrality_groups.pkl', 'rb'))
for key_c in centr_groups.keys():
    print(key_c, sorted(centr_groups[key_c], key=lambda i: i[1], reverse=True)[0:3])
info[info[3]=='10']
dcor_goals_continents_2 = pickle.load(open('distance_cor/goals/dcor_goals_continents_2.pkl', 'rb'))
dcor_goals_groups_2 = pickle.load(open('distance_cor/goals/dcor_goals_groups_2.pkl', 'rb'))
print(np.sum(np.sum(np.abs(dcor_goals_continents_2['World'] - dcor_goals_groups_2['Global North']))))
print(np.sum(np.sum(np.abs(dcor_goals_continents_2['World'] - dcor_goals_groups_2['Global South']))))
dict_cor_goals_continents_2 = pickle.load(open('distance_cor/goals/dict_cor_goals_continents_2.pkl', 'rb'))
dict_cor_goals_groups_2 = pickle.load(open('distance_cor/goals/dict_cor_goals_groups_2.pkl', 'rb'))
sorted(np.sqrt(dict_cor_goals_groups_2['Global South'].min_dcor.tolist()), reverse=True)<jupyter_output><empty_output><jupyter_text>Comparing Global South and Global North to World<jupyter_code># interlinkages
dict_cor_goals_continents_2 = pickle.load(open('distance_cor/goals/dict_cor_goals_continents_2.pkl', 'rb'))
dict_cor_goals_groups_2 = pickle.load(open('distance_cor/goals/dict_cor_goals_groups_2.pkl', 'rb'))
print(np.sum(np.abs(dict_cor_goals_groups_2['Global South'].min_dcor.to_numpy() - dict_cor_goals_continents_2['World'].min_dcor.to_numpy())))
print(np.sum(np.abs(dict_cor_goals_groups_2['Global North'].min_dcor.to_numpy() - dict_cor_goals_continents_2['World'].min_dcor.to_numpy())))
# centrality
centr_goals_continents = pickle.load(open('distance_cor/goals/centrality_continents.pkl', 'rb'))
centr_goals_groups = pickle.load(open('distance_cor/goals/centrality_groups.pkl', 'rb'))
centr_goals_groups['Global South']
goals_GS = []
central_GS = []
for centr in centr_goals_groups['Global South']:
    goals_GS.append(centr[0])
    central_GS.append(float(centr[1]))
    
goals_GN = []
central_GN = []
for centr in centr_goals_groups['Global North']:
    goals_GN.append(centr[0])
    central_GN.append(float(centr[1]))
    
goals_W = []
central_W = []
for centr in centr_goals_continents['World']:
    goals_W.append(centr[0])
    central_W.append(float(centr[1]))
print(np.sum(np.abs(np.asarray(central_GS) - np.asarray(central_W))))
print(np.sum(np.abs(np.asarray(central_GN) - np.asarray(central_W))))
DT.ODA.ALLD.CD
info[info[3]=='17']
dict_all['Kenya'].loc['DT.ODA.ALLD.CD']<jupyter_output><empty_output><jupyter_text>Pearson, Spearman and distance correlations<jupyter_code>from scipy.stats import pearsonr, spearmanr
noise = 0.1

A = np.random.uniform(-1,1,128).reshape(-1,1)
B = A**2 + np.random.uniform(-noise,noise,128).reshape(-1,1)
C = A**3 + np.random.uniform(-noise,noise,128).reshape(-1,1)
D = A**4 + np.random.uniform(-noise,noise,128).reshape(-1,1)
E = A + np.random.uniform(-noise,noise,128).reshape(-1,1)
fig, (ax1, ax2, ax3) = plt.subplots(1,3, figsize=(14,5))

ax1.set_title(r'$r$={:.2f} ($p$={:.2f}), $\rho$={:.2f} ($p$={:.2f}),'.format(pearsonr(A.reshape(-1),B.reshape(-1))[0], pearsonr(A.reshape(-1),B.reshape(-1))[1], spearmanr(A.reshape(-1),B.reshape(-1))[0], spearmanr(A.reshape(-1),B.reshape(-1))[1]) + '\n $dcor$={:.2f} ($p$={:.2f})'.format(dcor.distance_correlation(A,B), dcor.independence.distance_covariance_test(A,B, num_resamples=10000).p_value), fontsize=14)
ax1.set_xlabel('X', fontsize=14)
ax1.set_ylabel('Y', fontsize=14)
ax1.scatter(A, B);

ax2.set_title(r'$r$={:.2f} ($p$={:.2f}), $\rho$={:.2f} ($p$={:.2f}),'.format(pearsonr(A.reshape(-1),C.reshape(-1))[0], pearsonr(A.reshape(-1),C.reshape(-1))[1], spearmanr(A.reshape(-1),C.reshape(-1))[0], spearmanr(A.reshape(-1),C.reshape(-1))[1]) + '\n $dcor$={:.2f} ($p$={:.2f})'.format(dcor.distance_correlation(A,C), dcor.independence.distance_covariance_test(A,C, num_resamples=10000).p_value), fontsize=14)
ax2.set_xlabel('X', fontsize=14)
ax2.scatter(A, C);

ax3.set_title(r'$r$={:.2f} ($p$={:.2f}), $\rho$={:.2f} ($p$={:.2f}),'.format(pearsonr(A.reshape(-1),D.reshape(-1))[0], pearsonr(A.reshape(-1),D.reshape(-1))[1], spearmanr(A.reshape(-1),D.reshape(-1))[0], spearmanr(A.reshape(-1),D.reshape(-1))[1]) + '\n $dcor$={:.2f} ($p$={:.2f})'.format(dcor.distance_correlation(A,D), dcor.independence.distance_covariance_test(A,D, num_resamples=10000).p_value), fontsize=14)
ax3.set_xlabel('X', fontsize=14)
ax3.scatter(A, D);

plt.tight_layout()
plt.savefig('r_dcor.png')
print(r'$r$={:.2f} ($p$={:.2f}), $\rho$={:.2f} ($p$={:.2f}),'.format(pearsonr(A.reshape(-1),E.reshape(-1))[0], pearsonr(A.reshape(-1),E.reshape(-1))[1], spearmanr(A.reshape(-1),E.reshape(-1))[0], spearmanr(A.reshape(-1),E.reshape(-1))[1]) + '\n $dcor$={:.2f} ($p$={:.2f})'.format(dcor.distance_correlation(A,E), dcor.independence.distance_covariance_test(A,E, num_resamples=10000).p_value))<jupyter_output>$r$=0.99 ($p$=0.00), $\rho$=0.99 ($p$=0.00),
 $dcor$=0.99 ($p$=0.00)<jupyter_text>Differences between selected groupingsWe compute how different opposing socio-economic groupings are. We choose **Global South vs Global North** and **Low Income vs High Income** as the groupings of interest.<jupyter_code>cor_goals_groups_2['Global South'].to_csv('distance_cor/GS.csv')
cor_goals_groups_2['Global North'].to_csv('distance_cor/GN.csv')
cor_goals_groups_2['Low Income'].to_csv('distance_cor/LI.csv')
cor_goals_groups_2['Lower middle Income'].to_csv('distance_cor/LMI.csv')
cor_goals_groups_2['Upper middle Income'].to_csv('distance_cor/UMI.csv')
cor_goals_groups_2['High Income'].to_csv('distance_cor/HI.csv')
GS_GN = cor_goals_groups_2['Global South'] - cor_goals_groups_2['Global North']
LI_HI = cor_goals_groups_2['Low Income'] - cor_goals_groups_2['High Income']
LI_LMI = cor_goals_groups_2['Low Income'] - cor_goals_groups_2['Lower middle Income']
LI_UMI = cor_goals_groups_2['Low Income'] - cor_goals_groups_2['Upper middle Income']
LMI_UMI = cor_goals_groups_2['Lower middle Income'] - cor_goals_groups_2['Upper middle Income']
UMI_HI = cor_goals_groups_2['Upper middle Income'] - cor_goals_groups_2['High Income']

# saving files
GS_GN.to_csv('distance_cor/GS_GN.csv')
LI_HI.to_csv('distance_cor/LI_HI.csv')
LI_LMI.to_csv('distance_cor/LI_LMI.csv')
LI_UMI.to_csv('distance_cor/LI_UMI.csv')
LMI_UMI.to_csv('distance_cor/LMI_UMI.csv')
UMI_HI.to_csv('distance_cor/UMI_HI.csv')<jupyter_output><empty_output><jupyter_text>Which countries are closest to the average in each grouping?We want to find one country in each grouping that represents the grouping most accurately. To do so, we compute the average of each indicator over all countries in each grouping first, and calculate then the distance between each country in this grouping to these averages of the indicators. The distance between country $C$ and the average of the grouping $A$ can be defined as$$d(C, A) = \sqrt{\sum_{i=1}^n \sum_{t=1}^T (C_{it} - A_{it})^2},$$where $n$ is the number of indicators in this grouping, and $T$ the years. Here, the inner sum $\sum_{t=1}^T (C_{it} - A_{it})^2$ computes the squared Euclidean distance between any given indicator, before all $n$ indicators are summed up and sqaure-rooted.<jupyter_code>years = ['2000', '2001', '2002', '2003', '2004', '2005', '2006', '2007', '2008', '2009', '2010', '2011', '2012', '2013', '2014', '2015', '2016', '2017', '2018', '2019']
averages = {}
for group in ['Global North', 'Global South', 'Low Income', 'Lower middle Income', 'Upper middle Income', 'High Income']:
    print(group)
    #averages[group] = pd.DataFrame(index=indicators_values_i['Germany'].index, columns=years)
    averages[group] = pd.DataFrame(index=dict_all_std['Germany'].index, columns=years)
    
    for indicator in averages[group].index:
        for year in years:
            year_list = []
            for country in groups[group].dropna():
                #year_list.append(indicators_values_i[country].loc[indicator, year])
                year_list.append(dict_all_std[country].loc[indicator, year])
            averages[group].loc[indicator, year] = np.nanmean(year_list)
# saving averages
group_avg = open('distance_cor/goals/grouping_averages.pkl', 'wb')
pickle.dump(averages, group_avg)
group_avg.close()
# distance from each country in a grouping to the average of this grouping
distances = {}    
for group in ['Global North', 'Global South', 'Low Income', 'Lower middle Income', 'Upper middle Income', 'High Income']:
    print(group)
    distances[group] = {}
    
    for country in groups[group].dropna():
        e_list = []
        #for indicator in indicators_values_i[country].index:
        for indicator in dict_all_std[country].index:
            country_e = []
            averages_e = []
            j = 0    # counter

            for year in years:
                # we can only consider years with data available in both countries
                #if pd.isna(indicators_values_i[country].loc[indicator, year]) is False and pd.isna(averages[group].loc[indicator, year]) is False:
                if pd.isna(dict_all_std[country].loc[indicator, year]) is False and pd.isna(averages[group].loc[indicator, year]) is False:
                    #country_e.append(indicators_values_i[country].loc[indicator, year])
                    country_e.append(dict_all_std[country].loc[indicator, year])
                    averages_e.append(averages[group].loc[indicator, year])

                    j += 1

            # distance between country and grouping average per indicator
            if j > 0:
                e_list.append(distance.sqeuclidean(country_e, averages_e, w=1/j))
            else:
                e_list.append(np.nan)    # make NaN

        distances[group][country] = np.sqrt(np.nansum(e_list))
country_reps = {}
for group in ['Global North', 'Global South', 'Low Income', 'Lower middle Income', 'Upper middle Income', 'High Income']:
    country_reps[group] = min(distances[group], key=distances[group].get)
sorted_distances = {}
for group in ['Global North', 'Global South', 'Low Income', 'Lower middle Income', 'Upper middle Income', 'High Income']:
    sorted_distances[group] = sorted(distances[group].items(), key=lambda x: x[1], reverse=False)
    pd.DataFrame(sorted_distances[group]).to_csv('distance_cor/distances_average_{}.csv'.format(group))
for group in ['Global North', 'Global South', 'Low Income', 'Lower middle Income', 'Upper middle Income', 'High Income']:
    print(group)
    print('1st closest:', sorted_distances[group][0][1]/1000, sorted_distances[group][0][0])
    print('2nd closest:', sorted_distances[group][1][1]/1000, sorted_distances[group][1][0])
    print('3rd closest:', sorted_distances[group][2][1]/1000, sorted_distances[group][2][0])
    print('4th closest:', sorted_distances[group][3][1]/1000, sorted_distances[group][3][0])
    print('5th closest:', sorted_distances[group][4][1]/1000, sorted_distances[group][4][0])
    print('--------------------------')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17554444043383352, 'fraction_numerical': 0.03655958335481875, 'mean_word_length': 4.698893133509649, 'pattern_counts': {'":': 6, '<': 64, '<?xml version=': 0, '>': 79, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '186364', 'license': 'MIT', 'n_tokens_mistral': 24516, 'n_tokens_neox': 23083, 'n_words': 3886, 'path': '5_Additions.ipynb', 'repo_name': 'felix-laumann/SDG-dataset', 'size': '58184'}