Unsupervised learning on illicit trade data
Alice Lépissier
Bren School of Environmental Science and Management
Policy and media attention on illicit financial flows (IFF) has increased, with the recognition that Africa is a net creditor to the world.
 |
 |
|
| New York Times (2013) | Guardian (2015) | |
 |
 | |
| Guardian (2017) | Economist (2019) |
What is trade mis-invoicing?
- The deliberate mis-statement of price or quantity of internationally traded goods in invoices presented to customs
- Can occur at import or export
- Can result in an inflow or outflow of money
Motivations for trade mis-invoicing include:
- Evading tariffs
- Exploiting subsidy regimes
- Subverting forex and capital controls
- Hiding transfers of wealth
Mechanisms of mis-invoicing
Why does trade mis-invoicing matter?
- Outflows undermine the fiscal base and public spending
- Financing gap needed to meet the Sustainable Development Goals (SDGs)
- Combating trade mis-invoicing is crucial for the mobilization of domestic resources in the continent, and can catalyze sustainable development
 |
| Governance loop (credit: William Davis) |
Data source
- Trade mis-invoicing data-set of the United Nations Economic Commission for Africa
- Panel with $n=6,248,254$ of mis-invoiced trade between 179 reporting jurisdictions and 179 partner countries for years 2000-2016 and disaggregated commodities
- Citation: Lépissier, Alice, Davis, William, & Ibrahim, Gamal. (2019). Trade Mis-Invoicing Dataset (Version 1).
- The entire data-set
panel_results.Rdatawith ~6 million observations is available online. - The unit of observation is a reporter-partner-commodity-year tuple, where the data represent the illicit flow embedded in a transaction between a reporter $i$ and a partner $j$ for a commodity $c$ in year $t$.
- The full panel is ~2GB. This notebook works with smaller summary data-sets generated from the full panel. The code to generate these summary data-sets is available on the repository https://github.com/walice/Trade-IFF by running the script file
Compute IFF Estimates.R.
Methodology for calculating mis-invoiced trade
- We locate mis-invoicing in the discrepancies between reported trade flows and their mirrored statistics
- But not all observed discrepancies are due to illicit motives!
- Imports tend to be recorded on the basis of Cost of Insurance and Freight (CIF), while exports are recorded Free On Board (FOB)
Our approach
- Estimate the discrepancies between imports and mirror net exports as a function of both licit and illicit predictors
- Perform a harmonization procedure in order to generate a reconciled value that represents our best estimate of what the legitimate value of the trade should be on a FOB basis
- Calculate the IFF embedded in each transaction as the difference between the observed value and the reconciled value
- The panel can be aggregated over several dimensions, e.g. partner country, commodity, year. There are two methods to aggregate up the illicit flow:
- Net basis: simply add up all negative and positive values, so that inflows and outflows cancel each other out
- Gross Excluding Reversals (GER) basis: ignore all inflows, and sum up the positive values only across trading partners
Zoom in on Africa
- During 2000-2016, the continent lost on average \$83 billion a year in gross illicit outflows
- Net cumulative flows during that period were \$362 billion
Source: generated by Data Visualization.R in https://github.com/walice/Trade-IFF
Goals¶
This project will apply the following techniques to the data:
- Dimension reduction using Principal Components Analysis (PCA)
- Clustering
- Graph analysis
Data wrangling¶
- Mis-invoiced trade for countries (aggregated data)
- View 1: unit of observation = reporter-year; features = sectors
- View 2: unit of observation = sector-year; features = reporters
- Bilateral matrix of mis-invoiced trade (dyadic data)
- Unit of observation = reporter-partner-year
Mis-invoiced trade for countries by sectors
In [7]:
# Extract mis-invoicing in imports
IFF_Sector_Imp = IFF_Sector[['section', 'Imp_IFF_hi']]
IFF_Sector_Imp
Out[7]:
Mis-invoiced trade for dyads
In [12]:
# Extract mis-invoicing in imports
IFF_Dest_Imp = IFF_Dest[['partner.ISO', 'Imp_IFF_hi']]
IFF_Dest_Imp
Out[12]:
Metadata for countries
In [14]:
crosswalk = pd.read_excel("Data/crosswalk.xlsx").rename(columns={'Country': 'country'})
crosswalk.head()
Out[14]:
Auxiliary functions¶
In [15]:
def create_features(data, values, features, obs):
"""
Convert data-set in long format to wide and preserve information on year.
data: {IFF_Sector_Imp, IFF_Dest_Imp, IFF_Dest_AFR, ...}, as Pandas dataframe,
name of data-set from which to create feature space, must be in long format
values: {'Imp_IFF_hi', 'Exp_IFF_hi'}, as string, values that data-set will represent
features: {'reporter.ISO', 'section', 'partner.ISO'}, as string,
what to use as the feature space
obs: {'section', 'reporter.ISO'}, as string, what to use as the observation level
"""
features_data = data.pivot_table(values=values,
columns=features,
index=[obs, 'year'],
fill_value=0)
return features_data
In [17]:
def biplot_PCA(features_data, nPC=2, firstPC=1, secondPC=2, obs='reporter.ISO', show_loadings=False):
"""
Project the data in the 2-dimensional space spanned by 2 principal components
chosen by the user, along with a bi-plot of the top 3 loadings per PC, and color observations
by class label.
Args:
features_data: as Pandas dataframe, data-set of features
nPC: number of principal components
firstPC: integer denoting first principal component to plot in bi-plot
secondPC: integer denoting second principal component to plot in bi-plot
obs: string denoting index of class labels (in features_data)
show_loadings: Boolean indicating whether PCA loadings should be displayed
Returns:
plot (interactive)
pca_loadings (if show_loadings=True)
"""
# Run PCA (standardize data beforehand)
features_data_std = StandardScaler().fit_transform(features_data)
pca = PCA(n_components=nPC, random_state=234)
princ_comp = pca.fit_transform(features_data_std)
# Extract PCA loadings
cols = ['PC' + str(c+1) for c in np.arange(nPC)]
pca_loadings = pd.DataFrame(pca.components_.T,
columns=cols,
index=list(features_data.columns))
# Extract PCA scores
pca_scores = pd.DataFrame(princ_comp,
columns=cols)
pca_scores[obs] = features_data.reset_index()[obs].values.tolist()
pca_scores['year'] = features_data.reset_index()['year'].values.tolist()
score_PC1 = princ_comp[:,firstPC-1]
score_PC2 = princ_comp[:,secondPC-1]
# Generate plot data
if obs == 'reporter.ISO':
plot_data = pd.merge(pca_scores, obs_info, on=[obs, 'year'])
color_obs = 'reporter'
tooltip_obs = ['reporter', 'year', 'Income group (World Bank)', 'Country status (UN)']
else:
plot_data = pca_scores
color_obs = 'section'
tooltip_obs = ['section', 'year']
# Return chosen PCs to plot
PC1 = 'PC'+str(firstPC)
PC2 = 'PC'+str(secondPC)
# Extract top loadings (in absolute value)
# TO DO: use dict to iterate over
toploadings_PC1 = pca_loadings.apply(lambda x: abs(x)).sort_values(by=PC1).tail(3)[[PC1, PC2]]
toploadings_PC2 = pca_loadings.apply(lambda x: abs(x)).sort_values(by=PC2).tail(3)[[PC1, PC2]]
originsPC1 = pd.DataFrame({'index':toploadings_PC1.index.tolist(),
PC1: np.zeros(3),
PC2: np.zeros(3)})
originsPC2 = pd.DataFrame({'index':toploadings_PC2.index.tolist(),
PC1: np.zeros(3),
PC2: np.zeros(3)})
toploadings_PC1 = pd.concat([toploadings_PC1.reset_index(), originsPC1], axis=0)
toploadings_PC2 = pd.concat([toploadings_PC2.reset_index(), originsPC2], axis=0)
toploadings_PC1[PC1] = toploadings_PC1[PC1]*max(score_PC1)*1.5
toploadings_PC1[PC2] = toploadings_PC1[PC2]*max(score_PC2)*1.5
toploadings_PC2[PC1] = toploadings_PC2[PC1]*max(score_PC1)*1.5
toploadings_PC2[PC2] = toploadings_PC2[PC2]*max(score_PC2)*1.5
# Project top 3 loadings over the space spanned by 2 principal components
lines = alt.Chart().mark_line().encode()
for color, i, dataset in zip(['#440154FF', '#21908CFF'], [0,1], [toploadings_PC1, toploadings_PC2]):
lines[i] = alt.Chart(dataset).mark_line(color=color).encode(
x= PC1 +':Q',
y= PC2 +':Q',
detail='index'
).properties(
width=400,
height=400
)
# Add labels to the loadings
text=alt.Chart().mark_text().encode()
for color, i, dataset in zip(['#440154FF', '#21908CFF'], [0, 1], [toploadings_PC1[0:3], toploadings_PC2[0:3]]):
text[i] = alt.Chart(dataset).mark_text(
align='left',
baseline='bottom',
color=color
).encode(
x= PC1 +':Q',
y= PC2 +':Q',
text='index'
)
# Scatter plot colored by observation class label
points = alt.Chart(plot_data).mark_circle(size=60).encode(
x=alt.X(PC1, axis=alt.Axis(title='Principal Component ' + str(firstPC))),
y=alt.X(PC2, axis=alt.Axis(title='Principal Component ' + str(secondPC))),
color=alt.Color(color_obs, scale=alt.Scale(scheme='category20b'),
legend=alt.Legend(orient='right')),
tooltip=tooltip_obs
).interactive()
# Bind it all together
chart = (points + lines[0] + lines[1] + text[0] + text[1])
chart.display()
if show_loadings:
return pca_loadings
In [18]:
def scree_plot(features_data, show_explained_var=False):
"""
Create a cumulative scree splot and (optional) return the explained variance by each component.
features_data: as Pandas dataframe, the data-set on which to run PCA
show_explained_var: as Boolean, flag for whether to return explained variance
"""
features_data_std = StandardScaler().fit_transform(features_data)
pca = PCA(n_components=features_data_std.shape[1], random_state=234)
princ_comp = pca.fit_transform(features_data_std)
explained_var = pd.DataFrame({'PC': np.arange(1,features_data_std.shape[1]+1),
'var': pca.explained_variance_ratio_,
'cumvar': np.cumsum(pca.explained_variance_ratio_)})
# Adapted from https://altair-viz.github.io/gallery/multiline_tooltip.html
# Create a selection that chooses the nearest point & selects based on x-value
nearest = alt.selection(type='single', nearest=True, on='mouseover',
fields=['PC'], empty='none')
# The basic line
line = alt.Chart(explained_var).mark_line(interpolate='basis', color='#FDE725FF').encode(
alt.X('PC:Q',
scale=alt.Scale(domain=(1, len(explained_var))),
axis=alt.Axis(title='Principal Component')
),
alt.Y('cumvar:Q',
scale=alt.Scale(domain=(min(explained_var['cumvar']), 1)),
axis=alt.Axis(title='Cumulative Variance Explained')
),
)
# Transparent selectors across the chart. This is what tells us
# the x-value of the cursor
selectors = alt.Chart(explained_var).mark_point().encode(
x='PC:Q',
opacity=alt.value(0),
).add_selection(
nearest
)
# Draw points on the line, and highlight based on selection
points = line.mark_point().encode(
opacity=alt.condition(nearest, alt.value(1), alt.value(0))
)
# Draw text labels near the points, and highlight based on selection
text = line.mark_text(align='left', dx=5, dy=-5).encode(
text=alt.condition(nearest, 'cumvar:Q', alt.value(' '))
)
# Draw a rule at the location of the selection
rules = alt.Chart(explained_var).mark_rule(color='gray').encode(
x='PC:Q',
).transform_filter(
nearest
)
# Put the five layers into a chart and bind the data
out = alt.layer(
line, selectors, points, rules, text
).properties(
title='Cumulative scree plot',
width=500, height=300
)
out.display()
if show_explained_var:
return explained_var[['PC', 'var']]
PCA on feature space (for individual reporting countries)¶
Sector features¶
In [19]:
sector_features = create_features(IFF_Sector_Imp, 'Imp_IFF_hi',
features='section', obs='reporter.ISO')
sector_features
Out[19]:
Biplots¶
In [20]:
biplot_PCA(sector_features, 10, 1, 2, obs='reporter.ISO', show_loadings=True)
Out[20]:
Source: generated by Data Visualization.R in https://github.com/walice/Trade-IFF
In [22]:
biplot_PCA(sector_features, 10, 5, 6, obs='reporter.ISO')
Explained variance¶
In [23]:
scree_plot(sector_features, show_explained_var=True)
Out[23]:
Variance-stabilizing and normalizing transformations¶
In [24]:
# Plot distribution of illicit flow in each feature (i.e. sector)
fig, axes = joypy.joyplot(sector_features, colormap=plt.cm.viridis, figsize=(8,8),
title='Distribution of mis-invoicing across sectors');
In [28]:
sector_features_yeo = pd.DataFrame(sector_features_yeo,
index=sector_features.index,
columns=sector_features.columns)
fig, axes = joypy.joyplot(sector_features_yeo, colormap=plt.cm.viridis, figsize=(8,8),
title='Distribution of mis-invoicing across sectors (Yeo–Johnson transformation)');
In [30]:
biplot_PCA(sector_features_yeo, 10, 1, 2, obs='reporter.ISO')
Country features¶
In [32]:
country_features = create_features(IFF_Sector_Imp, 'Imp_IFF_hi',
features='reporter.ISO', obs='section')
country_features
Out[32]:
Biplots¶
In [33]:
biplot_PCA(country_features, 10, 1, 2, obs='section', show_loadings=True)
Out[33]:
Source: generated by Data Visualization.R in https://github.com/walice/Trade-IFF)
Source: generated by Data Visualization.R in https://github.com/walice/Trade-IFF)
Variance-stabilizing and normalizing transformations¶
In [44]:
biplot_PCA(country_features_log, 10, 1, 2, obs='section')
In [45]:
biplot_PCA(country_features_yeo, 10, 1, 2, obs='section')
PCA on bilateral trade matrix¶
In [46]:
partner_features = create_features(IFF_Dest_Imp, 'Imp_IFF_hi',
features='partner.ISO', obs='reporter.ISO')
partner_features
Out[46]:
Biplots¶
In [47]:
biplot_PCA(partner_features, 10, 1, 2, show_loadings=True)
Out[47]:
Source: generated by Data Visualization.R in https://github.com/walice/Trade-IFF
Intra-African illicit financial flows¶
In [61]:
partner_features_AFR = create_features(IFF_Dest_Imp_AFR, 'Imp_IFF_hi',
features='partner.ISO', obs='reporter.ISO')
partner_features_AFR
Out[61]:
In [65]:
biplot_PCA(partner_features_AFR, partner_features_AFR.shape[1], 1, 2)
Removing outliers¶
In [67]:
partner_features_AFR_noout = partner_features_AFR[~outlying]
In [68]:
biplot_PCA(partner_features_AFR_noout, 46, 1, 2)
Coloring by class label¶
In [72]:
biplot_PCA_classes(sector_features, 21, 1, 2, classes='Income group (World Bank)')
In [73]:
biplot_PCA_classes(partner_features, 46, 1, 2)