Censored Regression for Modelling International Small Arms Trading and its "Forensic" Use for Exploring Unreported Trades

TL;DR

This paper employs a censored regression model with network effects to analyze international small arms trade data, revealing endogenous network influences, path dependence, and enabling forensic analysis of unreported trades.

Contribution

It introduces a novel spatial autocorrelation model with endogenous network effects for small arms trade analysis, including a forensic application to detect unreported transactions.

Findings

Strong endogenous network effects identified

Evidence of path dependence in trade flows

Model effectively detects unreported trade tendencies

Abstract

In this paper we use a censored regression model to investigate data on the international trade of small arms and ammunition (SAA) provided by the Norwegian Initiative on Small Arms Transfers (NISAT). Taking a network based view on the transfers, we not only rely on exogenous covariates but also estimate endogenous network effects. We apply a spatial autocorrelation (SAR) model with multiple weight matrices. The likelihood is maximized employing the Monte Carlo Expectation Maximization (MCEM) algorithm. Our approach reveals strong and stable endogenous network effects. Furthermore, we find evidence for a substantial path dependence as well as a close connection between exports of civilian and military small arms. The model is then used in a "forensic" manner to analyse latent network structures and thereby to identify countries with higher or lower tendency to export or import than reflected in the data. The approach is also validated using a simulation study.

Figures

Tables

Table 1. Table 1 : Schematic representation of the evaluation scheme used in the simulation study. True conditions in the columns and Estimated in the rows. UR denotes under-reporting and UR ¯ ¯ UR \overline{\textbf{UR}} denotes censored observations that are not under-reported. Further abbreviations: True positive (TP), false positive (FP), false negative (FN) and true negative (TN).

TRUE
Estimated		UR	$\overline{\text{UR}}$	$\sum$
	UR	0	FP	FP
	$\overline{\text{UR}}$	0	TN	TN
	$\sum$	0	$0.75N$	$0.75N$

Table 2. (a) Classifier evaluation DGP1.

TRUE
Estimated		UR	$\overline{\text{UR}}$	$\sum$
	UR	0	FP	FP
	$\overline{\text{UR}}$	0	TN	TN
	$\sum$	0	$0.75N$	$0.75N$

Table 3. (b) Classifier evaluation DGP2.

TRUE
Estimated		UR	$\overline{\text{UR}}$	$\sum$
	UR	TP	FP	TP+FP
	$\overline{\text{UR}}$	FN	TN	FN+TN
	$\sum$	$0.1N$	$0.65N$	$0.75N$

Table 4. Table 2 : Different arms types included in the NISAT dataset with three digit arms category code, weapon type, subcategory and number of transfers in the dataset.

Code	PRIO Weapons Type	Subcategories
200	Small Arms
210		Pistols & Revolvers
230		Rifles/Shotguns (Military)
233		Assault Rifles
234		Carbines
235		Sniper Rifles
237		Semi-automatic Rifles (Military)
239		Shotguns (Military)
240		Machine Guns
243		Sub Machine Guns
245		Light Machine Guns
247		General Purpose Machine Guns
250		Military Weapons
260		Military Firearms
270		Machine Guns All Types
300	Light Weapons
310		Heavy Machine Guns <= 12.7mm
400	Ammunition
415		Small Arms Ammunition
417		Small Calibre Ammunition <= 12.7mm
418		Shotgun Cartridges
Source: nisat.prio.org.

Table 5. Table 3 : The 59 major exporting and importing countries of the small arms and ammunition dataset with ISO 3 country codes.

Country	ISO3 Code	Country	ISO3 Code	Country	ISO3 Code
Argentina	ARG	India	IND	Poland	POL
Australia	AUS	Indonesia	IDN	Portugal	PRT
Austria	AUT	Ireland	IRL	Romania	ROM
Belgium	BEL	Israel	ISR	Russia	RUS
Brazil	BRA	Italy	ITA	Saudi Arabia	SAU
Bulgaria	BGR	Japan	JPN	Serbia	SRB
Canada	CAN	Kenya	KEN	Singapore	SGP
Chile	CHL	South Korea	KOR	Slovenia	SVN
China	CHN	Kuwait	KWT	South Africa	ZAF
Colombia	COL	Latvia	LVA	Spain	ESP
Croatia	HRV	Lebanon	LBN	Sweden	SWE
Cyprus	CYP	Lithuania	LTU	Switzerland	CHE
Denmark	DNK	Malaysia	MYS	Thailand	THA
Egypt	EGY	Mexico	MEX	Turkey	TUR
Estonia	EST	Netherlands	NLD	Ukraine	UKR
Finland	FIN	New Zealand	NZL	Un. Arab Emirates	ARE
France	FRA	Norway	NOR	United Kingdom	GBR
Germany	DEU	Pakistan	PAK	United States	USA
Greece	GRC	Peru	PER	Uruguay	URY
Hungary	HUN	Philippines	PHL	-	-

Equations

$$\tilde{Y}=\sum _{k=1}^q{\rho }_k{W}_k\tilde{Y}+X\beta +ϵ,\text{ ϵ∼NN(0,σ2IN)}$$

$$\tilde{Y}=\bigg{(}\underbrace{I_{N}-\sum_{k=1}^{q}\rho_{k}W_{k}}_{\equiv A(\bm{\rho})}\bigg{)}^{-1}(X\beta+\epsilon)=\underbrace{(A(\bm{\rho}))^{-1}}_{\equiv B(\bm{\rho})}(X\beta+\epsilon)=B(\bm{\rho})(X\beta+\epsilon),$$

$$\begin{split}P(\tilde{Y}|X,\theta)&=\frac{1}{(2\pi\sigma^{2})^{\frac{N}{2}}}|A(\bm{\rho})|\exp\bigg{\{}-\frac{(A(\bm{\rho})\tilde{Y}-X\beta)^{\mathrm{T}}(A(\bm{\rho})\tilde{Y}-X\beta)}{2\sigma^{2}}\bigg{\}}.\end{split}$$

$$\tilde{Y}=\left(\begin{array}[]{c}\tilde{Y}_{o}\\ \tilde{Y}_{m}\\ \end{array}\right)\sim\mathcal{N}_{N}\left(\left(\begin{array}[]{c}\mu_{o}\\ \mu_{m}\\ \end{array}\right),\left(\begin{array}[]{cc}\Sigma_{oo}&\Sigma_{om}\\ \Sigma_{mo}&\Sigma_{mm}\\ \end{array}\right)\right),$$

$$B(\bm{\rho})X\beta=\left(\begin{array}[]{c}\mu_{o}\\ \mu_{m}\\ \end{array}\right)\text{, }B(\bm{\rho})(B(\bm{\rho}))^{\mathrm{T}}\sigma^{2}=\left(\begin{array}[]{cc}\Sigma_{oo}&\Sigma_{om}\\ \Sigma_{mo}&\Sigma_{mm}\\ \end{array}\right).$$

$$\footnotesize\begin{split}\ell_{\tilde{Y}_{m}|\tilde{Y}_{o}}(\theta)&=\log\bigg{(}\int_{(-\infty,c]^{N_{m}}}\frac{1}{\sqrt{(2\pi)^{N_{m}}|\Sigma_{m|o}|}}\exp\bigg{\{}-\frac{(U-\mu_{m|o})^{\mathrm{T}}\Sigma_{m|o}^{-1}(U-\mu_{m|o})}{2}\bigg{\}}dU\bigg{)},\end{split}\normalsize$$

$$\mathcal{M}=\{{\tilde{Y}}_m:{\tilde{Y}}_{m,1}<c,...,{\tilde{Y}}_{m,N_m}<c\}.$$

$${\tilde{Y}}^c\sim {\mathcal{N}}_{N_m}({\mu }_m+{\Sigma }_{mo}{\Sigma }_{oo}^{-1}({\tilde{Y}}_o-{\mu }_o),{\Sigma }_{mm}-{\Sigma }_{mo}{\Sigma }_{oo}^{-1}{\Sigma }_{om})$$

$$\tilde{Y}^{*}=\left(\begin{array}[]{c}\tilde{Y}_{o}\\ \mu^{c}_{m|o}\\ \end{array}\right)$$

$$\begin{split}S^{*}(\bm{\rho})&=\mathbb{E}_{\theta_{0}}[\tilde{Y}^{\mathrm{T}}(A(\bm{\rho}))^{\mathrm{T}}A(\bm{\rho})\tilde{Y}|\tilde{Y}_{o},X,\mathcal{M}]=\\ &=\text{tr}\bigg{(}(A_{mm}(\bm{\rho}))^{\mathrm{T}}A_{mm}(\bm{\rho})\Sigma_{m|o}^{c}\bigg{)}+(\tilde{Y}^{*})^{\mathrm{T}}(A(\bm{\rho}))^{\mathrm{T}}A(\bm{\rho})\tilde{Y}^{*}.\end{split}$$

$$\begin{array}{rl}Q(\theta |{\theta }_0)& =-\frac{N}{2}\log (2\pi {\sigma }^2)+\log (|A(\mathbf{\rho })|)-\frac{(S^{\ast }(\mathbf{\rho })-2{\beta }^{\mathrm{T}}{X}^{\mathrm{T}}A(\mathbf{\rho }){\tilde{Y}}^{\ast }+{\beta }^{\mathrm{T}}{X}^{\mathrm{T}}X\beta )}{2{\sigma }^2}.\end{array}$$

$$\begin{array}{rl}\hat{\beta }(\mathbf{\rho })& =(X^{\mathrm{T}}X{)}^{-1}{X}^{\mathrm{T}}A(\mathbf{\rho }){\tilde{Y}}^{\ast }\\ {\hat{\sigma }}^2(\mathbf{\rho })& =\frac{S^{\ast }(\mathbf{\rho })-{\tilde{Y}}^{\ast \mathrm{T}}(A(\mathbf{\rho }){)}^{\mathrm{T}}HA(\mathbf{\rho }){\tilde{Y}}^{\ast }}{N},\end{array}$$

$$\tilde{Q}(\bm{\rho}|\theta_{0})=\kappa+\log(|A(\bm{\rho})|)-\frac{N}{2}\log\bigg{(}S^{*}(\bm{\rho})-\tilde{Y}^{*\mathrm{T}}(A(\bm{\rho}))^{\mathrm{T}}HA(\bm{\rho})\tilde{Y}^{*}\bigg{)}.$$

$$\begin{array}{rl}\frac{\partial (A(\mathbf{\rho }){)}^{\mathrm{T}}A(\mathbf{\rho })}{\partial {\rho }_k}& =-W_k-W_k^{\mathrm{T}}+2{\rho }_k{W}_k^{\mathrm{T}}{W}_k+\sum _{l\ne k}{\rho }_l(W_k^{\mathrm{T}}{W}_l+W_l^{\mathrm{T}}{W}_k)=:R_k(\mathbf{\rho })\\ \frac{\partial (A(\mathbf{\rho }){)}^{\mathrm{T}}HA(\mathbf{\rho })}{\partial {\rho }_k}& =-H{W}_k-W_k^{\mathrm{T}}H+2{\rho }_k{W}_k^{\mathrm{T}}H{W}_k+\sum _{l\ne k}{\rho }_l(W_k^{\mathrm{T}}H{W}_l+W_l^{\mathrm{T}}H{W}_k)=:H_k(\mathbf{\rho }).\end{array}$$

$$R_k^{\ast }(\mathbf{\rho })=\text{tr}(R_{k,mm}(\mathbf{\rho }){\Sigma }_{m|o}^c)+({\tilde{Y}}^{\ast }{)}^{\mathrm{T}}{R}_k(\mathbf{\rho }){\tilde{Y}}^{\ast }$$

$$\frac{\partial \tilde{Q}(\mathbf{\rho }|{\theta }_0)}{\partial {\rho }_k}=-\text{tr}(B(\mathbf{\rho })W_k)-\frac{N}{2}\frac{R_k^{\ast }(\mathbf{\rho })-{\tilde{Y}}^{\ast \mathrm{T}}{H}_k(\mathbf{\rho }){\tilde{Y}}^{\ast }}{S^{\ast }(\mathbf{\rho })-{\tilde{Y}}^{\ast \mathrm{T}}(A(\mathbf{\rho }){)}^{\mathrm{T}}HA(\mathbf{\rho }){\tilde{Y}}^{\ast }}.$$

$$\begin{array}{rl}{\tilde{Y}}_{t,ij}& =X_{t,ij}^{\mathrm{T}}{\beta }_t+\underset{\text{Reciprocity}}{\underbrace{{\rho }_{t,1}{\tilde{Y}}_{t,ji}}}+\underset{\text{Exporter Effect}}{\underbrace{{\rho }_{t,2}\frac{1}{n-2}\sum _{u\ne j}{\tilde{Y}}_{t,iu}}}+\underset{\text{Importer Effect}}{\underbrace{{\rho }_{t,3}\frac{1}{n-2}\sum _{u\ne i}{\tilde{Y}}_{t,uj}}}+{ϵ}_{t,ij},\\ {ϵ}_{t,ij}& \sim \mathcal{N}(0,{\sigma }^2)\text{ for }i,j=1,...,n\text{, }i\ne j\text{, }t=1993,...,2014,n=59\text{ and }N=3422.\end{array}$$

$${\pi }_{t,ij}=\mathbb{P}(Z_{t,ij}>c_t|X_{t,ij},Z_{t,-ij};{\hat{\theta }}_t).$$

$${\Pi }_t=(I({\pi }_{t,ij}>J_t)).$$

$${\Gamma }_t=(I(Y_{t,ij}>0)).$$

$${\Omega }_t=({\omega }_{t,ij})={\Pi }_t-{\Gamma }_t$$

$$\begin{array}{rl}{\Omega }_t^{+}=& (I({\omega }_{t,ij}=1))\\ {\Omega }_t^{-}=& (I({\omega }_{t,ij}=-1)).\end{array}$$

$$\begin{array}{rl}\mathbf{\rho }& =(0.1,0.2,0.3{)}^{\mathrm{T}}\text{, }\beta =(1,2,3,4,5{)}^{\mathrm{T}}\text{, }p=5\text{, }n=20\text{, }N=380\\ X& \sim {\mathcal{N}}_p(1,I_p)\\ Z& \sim {\mathcal{N}}_N(B(\mathbf{\rho })X\beta ,B(\mathbf{\rho })B(\mathbf{\rho }{)}^{\mathrm{T}})\\ {\tilde{Z}}_{ij}& =I(Z_{ij}>q_{0.75}(Z))Z_{ij}\text{, for }i\ne j=1,...,n.\end{array}$$

$$\begin{array}{rl}{\Omega }_{\mathcal{T}}^{+}& =\sum _{t\in \mathcal{T}}{\Omega }_t^{+}\text{, }\\ {\Omega }_{\mathcal{T}}^{-}& =\sum _{t\in \mathcal{T}}{\Omega }_t^{-}.\end{array}$$

$${\ell }_{comp}(\theta )=-\frac{N}{2}\log (2\pi {\sigma }^2)+\log (|A(\mathbf{\rho })|)-\frac{(A(\mathbf{\rho })\tilde{Y}-X\beta {)}^{\mathrm{T}}(A(\mathbf{\rho })\tilde{Y}-X\beta )}{2{\sigma }^2},$$

$$\begin{array}{rl}\frac{\partial {\ell }_{comp}(\theta )}{\partial \beta }=& \frac{1}{{\sigma }^2}{X}^{\mathrm{T}}[A(\mathbf{\rho })\tilde{Y}-X\beta ]\\ \frac{\partial {\ell }_{comp}(\theta )}{\partial {\sigma }^2}=& -\frac{N}{2{\sigma }^2}+\frac{{\tilde{Y}}^{\mathrm{T}}(A(\mathbf{\rho }){)}^{\mathrm{T}}A(\mathbf{\rho })\tilde{Y}-2{\beta }^{\mathrm{T}}{X}^{\mathrm{T}}A(\mathbf{\rho })\tilde{Y}+{\beta }^{\mathrm{T}}{X}^{\mathrm{T}}X\beta }{2{\sigma }^4}\\ \frac{\partial {\ell }_{comp}(\theta )}{\partial {\rho }_k}=& -\text{tr}(B(\mathbf{\rho })W_k)\\ & -\frac{{\tilde{Y}}^{\mathrm{T}}[-W_k-W_k^{\mathrm{T}}+2{\rho }_k{W}_k^{\mathrm{T}}{W}_k+{\sum }_{l\ne k}{\rho }_l(W_k^{\mathrm{T}}{W}_l+W_l^{\mathrm{T}}{W}_k)]\tilde{Y}+2{\beta }^{\mathrm{T}}{X}^{\mathrm{T}}{W}_k\tilde{Y}}{2{\sigma }^2}.\end{array}$$

$$\begin{split}\frac{\partial^{2}\ell_{comp}(\theta)}{\partial\beta\partial\beta^{\mathrm{T}}}&=-\frac{1}{\sigma^{2}}X^{\mathrm{T}}X\\ \frac{\partial^{2}\ell_{comp}(\theta)}{\partial\beta\partial\sigma^{2}}&=-\frac{1}{\sigma^{4}}X^{\mathrm{T}}[A(\bm{\rho})\tilde{Y}-X\beta]\\ \frac{\partial^{2}\ell_{comp}(\theta)}{\partial\beta\partial\rho_{k}}&=-\frac{1}{\sigma^{2}}X^{\mathrm{T}}W_{k}^{\mathrm{T}}\tilde{Y}\\ \frac{\partial^{2}\ell_{comp}(\theta)}{\partial\sigma^{2}\partial\sigma^{2}}&=\frac{N}{2\sigma^{4}}-\frac{\tilde{Y}^{\mathrm{T}}(A(\bm{\rho}))^{\mathrm{T}}A(\bm{\rho})\tilde{Y}-2\beta^{\mathrm{T}}X^{\mathrm{T}}A(\bm{\rho})\tilde{Y}+\beta^{\mathrm{T}}X^{\mathrm{T}}X\beta}{\sigma^{6}}\\ \frac{\partial^{2}\ell_{comp}(\theta)}{\partial\rho_{k}\partial\sigma^{2}}&=\frac{\tilde{Y}^{\mathrm{T}}[-W_{k}-W_{k}^{\mathrm{T}}+2\rho_{k}W_{k}^{\mathrm{T}}W_{k}+\sum_{l\neq k}\rho_{l}(W_{k}^{\mathrm{T}}W_{l}+W_{l}^{\mathrm{T}}W_{k})]\tilde{Y}+2\beta^{\mathrm{T}}X^{\mathrm{T}}W_{k}\tilde{Y}}{2\sigma^{4}}\\ \frac{\partial^{2}\ell_{comp}(\theta)}{\partial\rho_{k}\partial\rho_{k}}&=-\text{tr}\bigg{(}B(\bm{\rho})W_{k}B(\bm{\rho})W_{k}\bigg{)}-\frac{\tilde{Y}^{\mathrm{T}}W_{k}^{\mathrm{T}}W_{k}\tilde{Y}}{\sigma^{2}}\\ \frac{\partial^{2}\ell_{comp}(\theta)}{\partial\rho_{k}\partial\rho_{l}}&=-\text{tr}\bigg{(}B(\bm{\rho})W_{l}B(\bm{\rho})W_{k}\bigg{)}-\frac{\tilde{Y}^{\mathrm{T}}(W_{k}^{\mathrm{T}}W_{l}+W_{l}^{\mathrm{T}}W_{k})\tilde{Y}}{2\sigma^{2}}.\end{split}$$

$$\frac{\partial \log (|A(\mathbf{\rho })|)}{\partial {\rho }_k}=-|A(\mathbf{\rho }){|}^{-1}\text{tr}[\text{adj}(A(\mathbf{\rho }))W_k]=-\text{tr}(B(\mathbf{\rho })W_k),$$

$$\frac{\partial\text{tr}(B(\bm{\rho})W_{k})}{\partial\rho_{l}}=\text{tr}\bigg{(}\frac{\partial B(\bm{\rho})}{\partial\rho_{l}}W_{k}\bigg{)}=-\text{tr}\bigg{(}B(\bm{\rho})\frac{\partial A(\bm{\rho})}{\partial\rho_{l}}B(\bm{\rho})W_{k}\bigg{)}=\text{tr}(B(\bm{\rho})W_{l}B(\bm{\rho})W_{k})$$

$$\frac{\partial \tilde{Q}(\mathbf{\rho }|{\theta }_0)}{\partial {\rho }_k}=-\text{tr}(B(\mathbf{\rho })W_k)-\frac{N}{2}\frac{R_k^{\ast }(\mathbf{\rho })-{\tilde{Y}}^{\ast \mathrm{T}}{H}_k(\mathbf{\rho }){\tilde{Y}}^{\ast }}{S^{\ast }(\mathbf{\rho })-{\tilde{Y}}^{\ast \mathrm{T}}(A(\mathbf{\rho }){)}^{\mathrm{T}}HA(\mathbf{\rho }){\tilde{Y}}^{\ast }}.$$

Full text

Keywords: Gravity Model; Latent Variable; Maximum Likelihood; Monte Carlo EM Algorithm; Network Analysis; Spatial Autocorrelation; Zero Inflated Data

1 Introduction

The Small Arms Survey Update 2018 indicates transfers of small arms in 2015 amounting to 5,7 billion (Holtom and Pavesi, 2018, p. 19) with a major share and highest increases in ammunitions (Holtom and Pavesi, 2018, p. 22). Given the often fatal consequences - civilian or military - of the availability of these arms for intrastate conflict and shootings as well as for interstate war, the absence of empirical evidence for supplier-recipient networks is surprising. A major reason behind this research gap are the notorious data deficiencies due to non-reporting and illicit trafficking (see Holtom and Pavesi, 2018, p. 29-46). Based on the only large-scale data base for small arms (Marsh and McDougal, 2016) we aim to analyse for the first time the small arms trading network. We integrate gravity models in a statistical network design to apply a forensic statistical analysis.

Starting with the seminal work of Tinbergen (1962), the gravity equation was quickly established as a valuable tool of empirical trade research. The success of the model stems from its intuitive interpretation as well as its surprisingly strong empirical validity, see e.g. Head and Mayer (2014). It is therefore not surprising that the concept was applied to all kinds of trade relations, including the international exchange of arms. An early example for these applications is the work of Bergstrand (1992). Although he doubted the suitability of the model for arms trade because of the strong political considerations in this area, the approach was taken up more recently. Akerman and Seim (2014) and Thurner et al. (2018) used the gravity model in order to explain whether Major Conventional Weapons (MCW) are exchanged. Martinez-Zarzoso and Johannsen (2017) rely on the framework of Helpman et al. (2008) to investigate the influence of economic and political variables on the so-called extensive and intensive margin of MCW trade. The interplay between oil imports and arms exports is determined using a gravity model in Bove et al. (2018). While the papers above focus on the exchange of MCW, in our paper we investigate transfers of small arms and ammunition (SAA) provided by the Norwegian Initiative on Small Arms Transfers (NISAT). This data is arguably even better suited for a gravity model since small arms are potentially less dependent on political decision making and many more trade occurrences are recorded.

We propose a network perspective on international SAA trade and conceptualize countries as nodes and transfers between them as directed, valued edges. Although gravity models are a standard tool for the analysis of dyadic data (Kolaczyk, 2009), endogenous network effects are rarely incorporated in these models. We do so by connecting the idea of gravity models with the spatial autoregressive (SAR) model adjusted to network data. Especially in sociology, SAR models are regularly used in a network context since the early eighties (Dow et al., 1982; Doreian et al., 1984; Doreian, 1989). They are called network autocorrelation models in this strand of literature. More recently, network autocorrelation models became popular in political science applications, see for example Franzese and Hays (2007), Hays et al. (2010) and Metz and Ingold (2017). Here, it is assumed that actors with certain characteristics are embedded in a network and this embedding leads to contagion and/or spillover effects transmitted through the edges that relate the actors (Leenders, 2002). Hence, one presumes that the characteristics of actors are correlated because a specific social, political or economic mechanism is connecting them. Note that the design of these models is different as compared to the usual set-up of gravity models since the outcome is related to the nodes, and the edges only represent indicators for node dependence. In this paper, we are interested in the dependencies among the transfers (instead of the actors) and account for outdegree, indegree, reciprocity and exogenous covariates. A similar model in a non-network context is the spatial gravity model (LeSage and Pace, 2008), that accounts for spatial dependence of the exporter, the importer as well as for the spatial importer-exporter dependencies.

Contrary to the typical structure of trade data we observe a high degree of reported non-trade in SAA. In other words, the trade network has a large percentage of zero entries. To accommodate the zero inflation problem we employ a censored SAR model that can be fitted using the Monte Carlo Expectation Maximization (MCEM) algorithm (Dempster et al., 1977; Wei and Tanner, 1990). There are already several similar EM-based approaches that have been pursued. For instance Suesse and Zammit-Mangion (2017) use the EM algorithm in spatial econometric models, Schumacher et al. (2017) apply an EM-based application to a censored regression model with autoregressive errors, and Vaida and Liu (2009) utilize EM estimation in a censored linear mixed effects model. In Augugliaro et al. (2018) a similar estimation procedure is used in the context of fitting a graphical LASSO to genetic networks.

While the model application per se provides new insights into SAA trading, the ultimate objective in this paper is to make use of the model to explore the validity of reported zero trades. This reflects a ”forensic” objective, i.e. we estimate, whether unreported trades are likely to have happened based on the fitted model. Despite this idea is in line with forensic statistics and forensic economics (Aitken and Taroni, 2004; Zitzewitz, 2012) our goal is apparently less ambitious. We do not aim to provide statistical evidence that some states are under-reporting but we do want to investigate potential under-reporting by utilizing the fitted network model.

This paper is organized as follows: after presenting the data in Section 2, we explain the model and show how to proceed with estimation and inference in Section 3. In Section 4 the results of the censored regression analysis are given and Section 5 provides the ”forensic” analysis, accompanied by a simulation study. Section 6 concludes the paper.

2 Data description

Since 2001, the Geneva-based Small Arms Survey specializes on documenting the international flows of the respective products. However, only the Norwegian Initiative on Small Arms Transfers (NISAT, see Marsh, 2017) provides truly relational data necessary for applying network analysis. The NISAT database contains relational information on the trade of small arms, light weapons and ammunition (see also Marsh and McDougal, 2016). This information is collected from different sources as described in Haug et al. (2002). Although NISAT represents the most reliable source of data regarding the exchange of small arms and light weapons, there is nevertheless an enormous amount of uncertainty inherent to arms trading data. This is especially true for light weapons where data quality and availability is partly very poor (Herron et al., 2011). Therefore, we restrict our analysis to small arms and the associated ammunition (SAA). See Table 2 in Annex 2 for the types of small arms and ammunition included in the dataset. Note, that the NISAT database also contains data on sporting guns, which we excluded from the dataset since we are particularly interested in the export of small arms with potential military value. Actually, we will rely on transferred sporting guns volumes later as a useful explanatory variable.

In the remaining dataset, more than 86 000 SAA transfers are recorded for the years 1992-2014, providing the exporting country, the importing country as well as the transferred arms category. The value of the export is measured in constant 2012 USD. In order to make estimation feasible, we restrict our analysis to a subnetwork and select those countries that account for the major share of the SAA trade activity. The resulting 59 countries (see Annex A.1, Table 3) account for $73\%-91\%$ (depending on the year) of the total transfer volume and have participated in arms trade at least once in each year under study. Hence, we investigate the ”core” of the international small arms trade network, balancing the trade-off between the number of countries included, the share of trade volume and the density of the subnetworks. In Figure 1, we show two binary networks for 1992 and 2014, with the countries represented as nodes and the arms transfers as directed edges among them.

In the left panel of Figure 2 we show the aggregated exports for the most important exporters United States (USA), Germany (DEU) and Italy (ITA) together with the exports of the 56 other countries (other). On the right hand side of Figure 2 we present the density, defined as the sum of existent edges divided by the number of potential edges. Although the network can be without doubt described as a dense one (as compared to the density of social networks), the density is smaller than $0.2$ in the beginning and remains below $0.4$ in the subsequent recent years.

3 Regression model

3.1 General model

Let $Y_{t}=(Y_{t,ij})\in\mathbb{R}^{n\times n}$ represent a network of transfers at the discrete time points $t=1,...,T$. At each time point $Y_{t}$ consists of $n$ nodes and $N=n(n-1)$ directed, continuous real valued edges, with diagonal elements $Y_{t,ii}$ left undefined. We set $\tilde{Y_{t}}=\text{vec}(Y_{t})\in\mathbb{R}^{N}$ as the row wise vectorization of $Y_{t}$, excluding the diagonal elements. In the following, we suppress the time index $t$ for ease of the notation and assume (after some suitable transformation) that in each time point $\tilde{Y}$ follows the autoregressive network model

[TABLE]

with $\beta$ being a $p$-dimensional parameter vector for the design matrix $X$. The matrices $W_{k}$ are row-normalized weight matrices representing linear endogenous network effects, with parameters $\rho_{k}$ as their strength. Model (1) is usually known as spatial autoregressive (SAR) model and we refer to LeSage and Pace (2009) for a more detailed discussion and to Lacombe (2004) or LeSage and Pace (2008) for similar models with multiple weight matrices. Standard software implementations that allow for a likelihood based estimation of the model are mostly restricted to the special case with $q=1$, for example in the R package spdep (Bivand et al., 2013; Bivand and Piras, 2015). The package tnam by Leifeld et al. (2017) allows for multiple weight matrices but is based on pseudo-likelihood estimation and therefore valid only if the weight matrices exclusively apply to exogenous covariates. Another possibility to estimate similar models is given by the package ARCensReg (Schumacher et al., 2017), initially designed to fit models with autoregressive errors. Because of the similar mathematical structure, the package could be used to fit models with spatially dependent errors known as Spatial Error Models (SEM). In the given case however, the network structure is assumed to influence the response directly which prevents us from using the package.

Model (1) can be rewritten as

[TABLE]

where the dependence on the $q$-dimensional parameter vector $\bm{\rho}=(\rho_{1},...,\rho_{q})^{\mathrm{T}}$ is made explicit for notational clarity. Similar as in Besag (1974) and given that all $N$ edges in the network are observed, their distribution is given by

[TABLE]

The parameter space of this model is restricted such that $A(\bm{\rho})$ is non-singular, which is ensured if the eigenvalues of $A(\bm{\rho})$ are real valued and greater than zero.

3.2 Censored regression model

The above model is not directly applicable to our data since a large proportion of the SAA trade values is zero, expressing no (reported) SAA trade between countries $i$ and $j$. We therefore adapt model (1) towards a utility model with censored observations. For each potential transfer from $i$ to $j$ there exists a utility of the transfer. This utility, however, only materializes in a transfer if it is higher than a certain threshold. Therefore, we assume that the probability model (2) applies to a network of partly observed latent variables, say $Z=(Z_{ij})$. The relation among $Y_{ij}$ and $Z_{ij}$ is given by $Y_{ij}=\max(c,Z_{ij})$ for $i,j=1,...,n$ and $c$ some threshold. Accordingly, we now set $\tilde{Y}=\text{vec}(Z)$ and label the observed utility $\tilde{Y}_{o}\in\mathbb{R}^{N_{o}}$ and the $N_{m}$ unobserved ones as $\tilde{Y}_{m}$. A reordering according to the observational pattern of $Y$ gives

[TABLE]

where $\tilde{Y}_{m}<c$ and $N=N_{o}+N_{m}$. Since the density of the network (see Figure 2) is roughly between $0.2$ and $0.4$ in all years, $N_{m}$ is always substantially larger than $N_{o}$. The mean-covariance structure is given by

[TABLE]

In the following, we will denote all reordered matrices in the notation with double subscripts, i.e. $A_{oo}$ refers to the submatrix of $A$ where only interactions of observed variables $\tilde{Y}_{o}$ enter.

3.3 Monte Carlo EM estimation

In order to estimate the unknown parameter vector $\theta=(\bm{\rho},\beta,\sigma^{2})\in\mathbb{R}^{q+p+1}$, we employ the EM algorithm (Dempster et al., 1977). The complete log-likelihood $\ell_{comp}(\theta)$ is simply derived from (2). We are interested in maximizing the observed log-likelihood $\ell_{obs}(\theta)=\ell_{\tilde{Y}_{o}}(\theta)+\ell_{\tilde{Y}_{m}|\tilde{Y}_{o}}(\theta)$, where the first part is simply the multivariate normal density of the observed transfers. The second part equals

[TABLE]

where $\mu_{m|o}$ and $\Sigma_{m|o}$ are the first and second conditional moments. Because $N_{m}$ is greater than $2\,000$ in each year, the observed log-likelihood is numerically hard to evaluate (and even more so to maximize) with state of the art software implementation. As a solution, we apply the EM algorithm and maximize $Q(\theta|\theta_{0})\coloneqq\mathbb{E}_{\theta_{0}}[\ell_{comp}(\theta)|\tilde{Y}_{o},X,\mathcal{M}]$ iteratively. The observation space is given by

[TABLE]

3.4 E-Step

The E-Step essentially boils down to calculating the first two moments of a multivariate normally distributed variable $\tilde{Y}^{c}$

[TABLE]

with restriction $\mathcal{M}$ from (3) applied to $\tilde{Y}^{c}$. Let those truncated moments be $\mu_{m|o}^{c}$ and $\Sigma_{m|o}^{c}$ and define

[TABLE]

as the vector that contains the observed values as well as the conditional expectation of the non-observed ones. Given the two moments, we can calculate the conditional expectation of the quadratic form (see Mathai and Provost, 1992):

[TABLE]

Then, the function to maximize in the M-step is given by

[TABLE]

In order to find the first and second moment of a truncated multivariate normally distributed variable, Vaida and Liu (2009) use the results of Tallis (1961) on the moment generating function to provide closed form expressions of the E-Step. This, however, is not practicable in our setting as (a) software implementations of a multivariate normal distribution function are overstrained by the high dimension of our problem (the standard package in R, mvtnorm by Genz et al. (2016) is not able to process dimensions higher than $1\,000$) and (b) as noted by Schumacher et al. (2017), even if the distribution function could be evaluated, the closed form solution is computationally very expensive which leads to infeasible convergence times in applications with a high number of non observed values. The same is true for the direct calculation using the moment generating function implemented in R by Wilhelm et al. (2012).

A practicable alternative consists in using the Monte Carlo EM (MCEM) algorithm (Wei and Tanner, 1990) where intractable expectations are replaced by sample based approximations. In our specific case we use the R package TruncatedNormal by Botev (2017) in order to draw from the truncated multivariate normal distribution. An alternative would be to enrich the E-Step with a stochastic approximation step (SAEM algorithm, see Schumacher et al., 2017 for a detailed description) which reduces the number of simulations needed and is very efficient if the M-step is faster then the E-Step. In our specific application, the computational bottleneck comes with the M-Step and simulations showed that the SAEM converges more slowly than the MCEM algorithm.

3.4.1 M-Step

It is numerically more efficient to reduce the log-likelihood to a profile log-likelihood by first maximizing with respect to $\beta$ and $\sigma^{2}$ and then with respect to $\bm{\rho}$. Using the derivatives of (7) with respect to $\beta$ and $\sigma^{2}$ and defining $\hat{\beta}(\bm{\rho})$ and $\hat{\sigma}^{2}(\bm{\rho})$ as the solutions of the score equations as functions of $\bm{\rho}$ it follows that

[TABLE]

where $H=X(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}$ is the hat matrix. With $\kappa$ being a constant we can write the profiled function $\tilde{Q}(\cdot)$ as

[TABLE]

The expressions $(A(\bm{\rho}))^{\mathrm{T}}A(\bm{\rho})$ and $(A(\bm{\rho}))^{\mathrm{T}}HA(\bm{\rho})$ have derivatives

[TABLE]

Now define

[TABLE]

which gives

[TABLE]

Iteration between the E- and the M-step provide the final estimate $\hat{\theta}$. The variance of $\hat{\theta}$ can be calculated using Louis (1982) formula with more details on the practical implementation provided in the Supplementary Material.

4 Application to the data

4.1 Covariates

Considering model (1) we need to specify the two major components of the model, namely (a) the covariates included in matrix $X$ and (b) the network related correlation structure.

Node Specific Variables: Following standard applications (Ward et al., 2013; Head and Mayer, 2014; Egger and Staub, 2016; Thurner et al., 2018), we control for the logarithmic real GDP in constant 2010 USD as a measure for the market size of the exporting and importing country. The data are provided by the World Bank (2017). For the two years 1993-1994 no reliable GDP data are available for Serbia, Croatia, Estonia, Latvia, Lithuania and Slovenia, we therefore assume that the GDP remained constant in the first three years for this countries. In order to control for the potential influence of intrastate conflicts we insert a binary variable that is one if there is an intrastate conflict in the receiving country in the respective year and zero otherwise. The corresponding data is available from the webpage of the Uppsala Conflict Data Program (UCDP, 2019).

Edge Specific Distance Measures: Because of the strong empirical evidence that geographic distance is a relevant factor in trade (Disdier and Head, 2008), we control for the logarithmic distance between capital cities in kilometres (Gleditsch, 2013). In recent applications of the gravity model to arms trade (Akerman and Seim, 2014; Martinez-Zarzoso and Johannsen, 2017; Bove et al., 2018; Thurner et al., 2018) it is argued that political distance measures in terms of regime dissimilarity must also be inserted in the gravity equation. We use the absolute difference of the polity IV index (Marshall, 2017) between two countries, ranging from 20 (highest ideological distance) to zero (no ideological distance). Additionally, we include a dummy variable for formal alliances between the exporting and importing country, being one if the two countries have a formal alliance. The data is available from Correlates of War Project (2017) until 2012 and we assume that the alliances stay constant for the years 2013 and 2014.

Edge Specific Trade Measures: We control for lagged logarithmic SAA transfers by smoothing the past observed trade volume using a five-year moving average. In the years with less then five lagged periods available, the moving average is shortened accordingly. We call this path dependency, leading to inertias that arises because of diminishing transaction costs, trust relations, security aspects and potentially interoperability, and is a very important determinant in the MCW trade network (Thurner et al., 2018).

Additionally, we enrich the model with a five year moving average of logarithmic civilian weapon transfers. The intuition behind that is that exports of SAA for military usage and civilian usage might be correlated. This is plausible because countries that export massive amounts of civilian arms also have the capabilities to produce military arms. The data is also provided by NISAT (Marsh and McDougal, 2016). Furthermore it seems plausible that there is a connection between the volume of small arms traded and the volume of MCW. MCW transfers are recorded by the Stockholm International Peace Research Institute (SIPRI) and measured in so called trend indicator values (TIV). This measure represents the military value and the production costs of the transferred products. For detailed explanation of the data and the TIV see SIPRI (2017b, a) and Holtom et al. (2012). We use a dummy variable that is one if there was an MCW transfer from country $i$ to $j$ in the actual year or in the four preceding years, zero otherwise. Additionally, we use the logarithmic sum of the exported TIV volumes in the actual year and the four preceding ones.

4.2 Network structure

Next we specify the network specific effects represented by matrices $W_{k}$ in model (1). We include three effects which are explained subsequently and visualized in Figure 3.

Reciprocity: The reciprocity effect measures whether the export volume from country $i$ to country $j$ increases in the export volume from $j$ to $i$. In the given context it is a plausible assumption that countries tend to specialize in certain types of small arms and/or ammunition and therefore complement each other with their products. Mutual trade is likely to be encouraged by political partnerships and indicates strategic elements, induced by bilateral agreements. The measure is also investigated in the context of commercial trade (e.g. Garlaschelli and Loffredo, 2005; Barigozzi et al., 2010;Ward et al., 2013). In the arms trade literature, reciprocity is specified by Thurner et al. (2018), with the finding that this is rather unusual in the context of MCW.

Exporter and Importer Effect: The exporter and the importer effect have their analogies in binary networks and can be interpreted as the valued versions of the outdegree and the indegree. The coefficient of the exporter effect measures whether the transfers going out from a certain exporter $i$ are correlated. A positive effect indicates the presence of ”super-exporters”. Contrary, the importer effect measures whether the imports of a certain importer $j$ are related, with a positive effect indicating ”super-importers”. The degree structure is a crucial feature of the SAA network because a rather small number of countries accounts for the major share of the trade volume, while a small share of (potentially identical) importing countries accounts for a great amount of the import volume.

Before fitting the model we apply the natural logarithm to the data. This is necessary, because in its raw form the data is strongly skewed with a long tail. In the Supplementary Material the distribution of the log-transformed response is investigated

Hence, if the original trade matrices are given by $Y_{t}=(Y_{t,ij})$, the elements of $\tilde{Y}_{t}$ are given by $\tilde{Y}_{t,ij}=\log(Y_{t,ij})$ if $Y_{t,ij}>0$ and are not defined if $Y_{t,ij}=0$. Furthermore, we define $d_{t}=\min(\{Y_{t,ij}>0\})$ as the lowest strictly positive value in the network at year $t$ and set $c_{t}=\log(d_{t})$. That is, the threshold $c_{t}$ is defined such that at a given time point $t$ all transfers below the smallest observed log-transformed transfer in that sample are censored. Utility below the threshold $c_{t}$ implies that no transfer was carried out or was not recorded. Furthermore, we allow for time-varying coefficients by estimating each time-period separately. This relaxes the unrealistic assumption of time-constant effects for more than $20$ years and reduces the computational effort. Given these specifications, the final model is now given by

[TABLE]

4.3 Results: Coefficients

In Figure 4 the time series of coefficients are plotted against time for the years 1993-2014. The shaded areas around the coefficients give two standard error bounds and the colouring of the point estimates reflects the respective significance level, the zero-line is depicted by solid black.

Covariates: The exogenous covariates in the first row and in the second row on the right represent the standard gravity variables logarithmic GDP of the exporter and the importer as well as the geographic distance between them (second row, right panel). Overall, the expected results of the gravity equation hold, except for the logarithmic GDP of the exporter that tends to be insignificant and is close to zero in the most recent years. This is an interesting result, because it highlights the fact, that market size is not a prerequisite for producing and exporting internationally competitive SAA. This finding is in stark contrast to the insights on MCW by Thurner et al. (2018).

The strong negative effect of the geographic distance is significant in all years. Again, this is different as compared to MCW transfers where geopolitical strategy disregards distance.

Regarding the political security measures we find that the presence of a conflict in the importing country (second row, left panel) has a mostly positive but seldom significant effect while the coefficient on the dissimilarity of political regimes (third row, left panel) is mostly negative but also often insignificant. The coefficient on the dummy variable for formal alliances (third row, right panel) is positive in the beginning but almost permanently insignificant from 2001 on.

The large and consistently significant coefficients of the lagged moving average (fourth row, left panel) illustrates an important feature of the network, namely path inertia. Intensive transfer relationships in the past, strongly increase the export volume in the present. Similarly we find a strong connection between exporting civilian and military arms (fourth row, right panel). Looking at the relation between SAA and MCW trade we find that having traded MCW (fifth row, left panel) in the actual year or the in the four preceding ones has a strong positive effect - at least until the last two years. However, at the same time the effect of the logarithmic sum of the TIV values (fifth row, right panel) has a negative and mostly significant effect. I.e. rather small transfers of MCW tend to coincide with SAA exports while dyads with high amounts of MCW exchange tend to transfer small arms to a relatively lower degree.

Network Structure: On the right panel in the sixth row of Figure 4 the coefficients for reciprocity are shown. The coefficients remain almost constant and positive with values between $0.02$ and $0.06$. As the coefficients often changes from significant to insignificant we infer that there is at least a tendency that mutuality increases the volume of arms exchanged.

The strongest endogenous effect is the exporter effect (bottom row, left panel) with coefficients that are consistently positive and significant. This indicates that the transfers stemming from the same exporter are indeed highly correlated and reflects the existence of ”super-sellers” like the United States, Germany, Brazil or Italy. On the other hand, we also find a stable, positive and significant importer effect (bottom row, right panel). The fact that the two coefficients on the exporter and the importer effect are much higher than the reciprocity effect provides structural information about heterogeneity in the network. Being a strong exporter or sending to a strong importer increases the export volume more than simply having imported high amounts from the respective partner.

5 ”Forensic” statistical analysis

5.1 Under- and over-reporting

Our model rests on the assumption that the SAA network is determined by a latent utility network $Z_{t}$. Based on the joint distribution (2) we can in fact estimate the probability of $Z_{t,ij}$ being greater than the threshold $c_{t}$, given the covariates, the endogenous effects and the rest of the network. In order to do so, let $Z_{t,-ij}$ represent the $(N-1)$-dimensional vector that contains the realized and the expected values of the latent variables, except the entry that corresponds to the transfer from $i$ to $j$. Because we are interested whether some latent transfers could have realized according to the model, we form the expectations without the restriction that the latent transfers must be smaller than $c_{t}$. Based on this, we define the conditional probability of a specific latent transfer being greater than the threshold by

[TABLE]

By construction (see the Supplementary Material for the derivation), $\pi_{t,ij}$ is high for transfers that are observed in the dataset ($Y_{t,ij}>0$) and small for transfers that are not observed ($Y_{t,ij}=0$). However, we may calculate a high value of $\pi_{t,ij}$, i.e. a high probability for a realized transfer of arms, despite the data actually indicates $Y_{t,ij}=0$. We propose to consider this as potential under-reporting. Such a zero-record can happen due to random fluctuation, factors beyond the model as for example historical relationships, or because de-facto existent transfers have not been reported. Vice versa, we may obtain a low value of $\pi_{t,ij}$ although $Y_{t,ij}>0$. We label this as over-reporting. This label is not intended to suggest that potentially over-reported transfers in fact never happened, but highlights transfers where our model attaches a lower level of latent utility than manifested in the data. Naturally, our main ”forensic” interest is in uncovering potential under-reporting.

Apparently, this requires the fixation of a threshold value for the probabilities. Based on Receiver-Operating Characteristic (ROC) curves, an optimal threshold value $J_{t}$ can be found using Youden’s $J$ statistic (Youden, 1950). This value is optimal in the sense that it allows for a separation such that both, sensitivity and specificity are maximized. This defines the binary network

[TABLE]

This network is now set into relation with the observed binary SAA trade

[TABLE]

Comparing $\Pi_{t}$ and $\Gamma_{t}$, we can define the ”forensic” network

[TABLE]

which in turn creates two new binary networks

[TABLE]

For $\omega_{t,ij}=1$, the model predicted a transfer that is not present in the dataset, and for $\omega_{t,ij}=-1$, the model did not predict an actual transfer. Following our convention from above we label $\Omega_{t}^{+}$ as the under-reporting network and to $\Omega_{t}^{-}$ as the over-reporting network of unpredicted but realized transfers.

5.2 Simulation study of ”forensic” power

Before we apply our model in a ”forensic” matter to identify transfers with potential under-reporting we demonstrate the behavior of the model in a simulation study to explore its detection properties. We use two different settings in order to investigating how well the proposed approach identifies under-reporting. The first setting builds on the following Data Generating Process (DGP1)

[TABLE]

Here, $q_{0.75}(Z)$ denotes the 75% quantile and we are censoring the network towards an observed density of $0.25$. Note, that DGP1 is not subject to under-reporting and all censored responses are in fact below the censoring threshold. The results of running DGP1 $100$ times and applying the estimation procedure are summarized in the Supplementary Material, indicating that the expected values approximate the latent variables very well and that we are able to find unbiased estimates despite the enormous amount of censoring.

In order to validate the forensic power of the model, we run a second experiment (DGP2), being a modified version of DGP1. To be precise, we are censoring again 75% of the observations but only 65% correspond to the lowest ones, while the remaining 10% are randomly selected among the observations that are in fact higher than the threshold $q_{0.65}(Y)$. This share of observations represents the under-reporting. Again we run DGP2 100 times.

In order to make the following evaluation transparent, we represent the evaluation scheme (e.g. Fawcett, 2006) for both DGPs in Table 1. On the left hand side, we regard the simulation without under-reporting (DGP1). In this setting we can investigate the false positive rate (FPR), being the sum of the false positives (FP) relative to the number of all observations that are in fact not under-reported. A low value for this measures means in DGP1 that in a setting without under-reporting a low share is classified as under-reporting. In DGP2, the measure tells us whether including true under-reporting in the simulation leads to an increase of misclassified under-reporting. The corresponding results are visualized on the top panel of Figure 5. The FPR shows a higher variability in DGP2 (right panel) and is slightly higher as compared to DGP1 (left panel). However, the results provide evidence for an overall low FPR in both setting.

Furthermore, DGP2 allows to evaluate the share of under-reported observations that is identified. This is assessed based on the true positive rate (TPR) and shown on the bottom left panel of Figure 5. In fifty percent of all simulation runs we are able to identify at least 95% of the falsely censored observations and even in the simulation runs with the worst performance, the TPR are does not fall below 74%. Additionally, we investigate the False discovery rate (FDR) that relates the observations that are wrongly classified to be under-reporting to the sum of all observations that are classified for under-reporting. A low value for this measure provides evidence, whether the model is able to keep the number of potential over-reporting that are in fact not under-reporting low. The corresponding results are shown in the south-east panel of Figure 5. We find a median share of less than 26% to be classified incorrectly.

5.3 ”Forensic” analysis of arms trade data

We now turn back to the data and provide the development of the densities for the latent networks in Figure 6. In the real data, over- and under-reporting is certainly not random but potentially clustered among countries. We therefore evaluate node (i.e. country) specific network topologies of $\Omega_{t}^{+}$ and $\Omega_{t}^{-}$ for each year and summarize the information in box-plots for each country, ordered according to the median of the respective feature. This is shown in Figure 7 for potential under-reporting and in Figure 8 for over-reporting. In the first row, we represent the Eigenvector centrality scores. This measure is undirected and constructed such that the centrality of each country is proportional to the sum of the centralities of its trading partners. Hence, countries with high scores have many potentially under-reported (over-reported) import- and export-relations with many other countries that themselves have many under-reported (over-reported) import- and export-relations, see e.g. Csardi and Nepusz (2006). In the middle row, we present the outdegree, that is the number of potentially under-reported (over-reported) exports for a country. The bottom row in Figures 7 and 8 gives the indegree, that is the number of potentially under-reported (over-reported) imports. All measures are scaled to take values between [math] and $1$. Countries at the right hand side in the plots of Figure 7 are potentially under-reporting and in Figure 8, the right hand side of the plots mirrors high over-reporting.

To detect persistent patterns in the networks on a dyadic level, we check whether potential under-reporting or over-reporting occurs frequently, i.e. counting instances of $\omega_{t,ij}^{+}=1$ and $\omega_{t,ij}^{-}=1$ for $t\in\mathcal{T}=\{1993,...,2014\}$. Denote the aggregated ”forensic” networks as

[TABLE]

We look at the distribution of elements of $\Omega_{\mathcal{T}}^{+}$ and $\Omega_{\mathcal{T}}^{-}$, which is plotted in Figure 9. On the horizontal axis we show the possible values of the matrix entries, that is the number years where transfers in the forensic networks occur. This ranges from 1 (potential under-reporting or over-reporting in one year) to 22 (potential under-reporting or over-reporting in all years). The maximum entry of $\Omega_{\mathcal{T}}^{+}$ is thereby less then 22, namely 21, while the maximum value of $\Omega_{\mathcal{T}}^{-}$ is 15. On the vertical axis of Figure 9 we show the frequency of the entries of $\Omega_{\mathcal{T}}^{+}$ and $\Omega_{\mathcal{T}}^{-}$. Apparently for ”forensic” purposes, large values of $\Omega_{\mathcal{T}}^{+}$ are of particular interest, since they report pairs of countries which are likely to under-reporting.

The line in solid black with the ”+” symbols represents $\Omega^{+}_{\mathcal{T}}$ and the line in grey with the ”-” symbols represents the under-reporting network. Additionally, we indicate for both networks the pairs of countries (i.e. sender and receiver) which are of particular interest for ”forensic” purposes. This means for example for an element of $\Omega_{\mathcal{T}}^{+}$ that has value 21, that the respective transfer from $i$ to $j$ is one of the four transfers appeared that appeared $21$ times in the under-reporting network.

Under-reporting networks $\Omega_{t}^{+}$: Looking at the Eigenvector centrality scores of Figure 7 on the top provides conclusive results about countries that are central in the network series $\Omega_{t}^{+}$. Among the countries where arms transfers are potentially under-reported, we find many Western European countries such as Belgium (BEL), Sweden (SWE), France (FRA), Spain (ESP) and Denmark (DNK). However, the list of presumed under-reporting is headed by Russia (RUS) and Turkey (TUR) but also Brazil (BRA), Israel (ISR) and China (CHN) have high scores. These countries also play a dominant role in Figure 9. In particular, exports from Brazil (BRA) to Russia (RUS), Hungary (HUN), Ukraine (UKR), China (CHN) and Japan (JPN) are likely to be frequently under-reported. Similarly, exports from Russia (RUS) to Cyprus (CYP), Denmark (DNK), Ireland (IRL), Turkey (TUR), Sweden (SWE), Portugal (PRT) and Greece (GRC) are listed. We also find imports of Israel (ISR) from Finland (FIN), Belgium (BEL) and Austria (AUT) as well as exports of Belgium (BEL) to India (IND), Ukraine (UKR), Russia (RUS), Lebanon (LEB), Colombia (COL) and China (CHN).

Over-reporting networks $\Omega_{t}^{-}$: Among the twelve countries with the highest Eigenvector centrality in Figure 8 is Croatia (HRV) as the only European country. There are however many countries from Asia such as Singapore (SGP), India (IND), Thailand (THA), Indonesia (IDN), Malaysia (MYS), South Korea (KOR) and Philippines (PHL). Furthermore, South Africa (ZAF), New Zealand (NZL), Australia (AUS) and Israel (ISR) are among the countries where trade activity is often over-reported. For Asian countries as well as for Australia and New Zealand this might mirror the fact that those countries export many SAA to Europe and the United States despite the strongly negative distance effect of the model. Furthermore, this network is very likely to be driven by bilateral agreements and historical developments not covered by the covariates. See for example in Figure 9 the number of over-reporting related to the Baltic countries Estonia (EST) Lithuania (LTU) and Latvia (LVA).

It remains to be emphasized that the constructions of the forensic networks relies on our model with corresponding assumptions and admittedly high degrees of uncertainty. As a consequence, it does not allow for definite statements about actual hidden transfers. However, many of the dyads listed in Figure 9 indeed have either traded massive amounts of civilian arms (e.g. AUT-CHN, BRA-RUS, RUS-BEL, RUS-DNK) or had frequent MCW trade relations (e.g. RUS-CYP) but almost no documented small arms transfers for military usage. Additionally, many of the countries that take central positions in the forensic networks are known for not being very transparent with respect to their SAA exports and imports (see e.g. the small arms transparency barometer).

6 Conclusion

In this paper we have modelled the volumes of international transfers of small arms and ammunition for the years 1992-2014 based on data provided by NISAT. As an analytical tool we combined the gravity model of trade with a modified SAR model that allows to enrich the analysis by endogenous network dependencies, accounting for exporter-related, importer-related and reciprocal dependency among the transfers in the network. Using a censored normal regression model we are able to include information provided by zero-valued transfers. The infeasible likelihood of the censored model is maximized using a Monte Carlo EM algorithm. The fitted model shows strong and stable endogenous network effect, especially related to the sender effect and the receiver effect but also some evidence for reciprocity. Additionally, we find a high coefficient on path dependency and a close connection to the exports of civilian small arms. Conditional on that, the classical gravity hypothesis is confirmed with respect to the GDP of the importer and physical distance but only exceptionally with respect to political distance measures and the GDP of the exporter. This contrasts with the MCW network where distance plays no role, where political similarity and GDP of the exporter have a strong impact (see Thurner et al., 2018). Actually, this difference is plausible, as the technological requirements for the production of small and ammunition are relatively low, and strategic considerations of world-wide acting countries make geographic distances a negligible factor for MCW trade.

Building on our latent utility framework we were able to explore latent utility networks. With the construction of under-reporting and over-reporting networks we perform for the first time a forensic approach in this area highlighting especially potentially under-reported exports of Russia and Turkey. We refrain, of course, from making too far-reaching assertions. Note that we do not claim to provide unambiguous claims for intentional false reporting. However, we demonstrate that some zero entries in the SAA trading network tend to be not plausible.

Acknowledgements

We would like to thank Nic Marsh from NISAT for his useful comments and explanations on the data.

Appendix A Annex

A.1 Descriptives

Appendix B Supplementary Material

B.1 Derivatives of the complete log-likelihood

The complete log-likelihood is given by

[TABLE]

with score vector

[TABLE]

And the corresponding Hessian results in

[TABLE]

Where we use Jacobi’s formula (see Magnus and Neudecker, 1988) that allows to express the derivative of a matrix determinant in terms of the derivative of the matrix and its adjugate ($\text{adj}(\cdot)$). Resulting in

[TABLE]

for the third equation in (12). The differentiation of the trace

[TABLE]

is used for the sixth and seventh equation in (13).

B.2 Practical Implementation of the Algorithm

The gradient

[TABLE]

can be used to maximize

[TABLE]

by applying the BFGS optimization routine (see Broyden, 1970, Fletcher, 1970, Goldfarb, 1970 and Shanno, 1970). The implementation of the BFGS algorithm in R (R Core Team, 2016) is provided by the base function optim. More computational stability for the maximization of equation (15) is reached by defining $\lambda=(\lambda_{1},...,\lambda_{N})^{\mathrm{T}}$ as the vector of eigenvalues of $A(\bm{\rho})$ and replacing $\log(|A(\bm{\rho})|)$ by $\sum_{r=1}^{N}\log(\lambda_{r}(\bm{\rho}))$ in equation (15), see Bivand and Piras (2015). The starting value for the algorithm can be found by using a maximum pseudolikelihood estimate (MPLE), using $W_{1}Y,...,W_{q}Y$ as exogenous covariates in a censored regression model, provided by the R package censReg (Henningsen, 2013). Since the observed log-likelihood cannot be evaluated, we define $\hat{\theta}$ as the solution of the maximization problem if $(\hat{\theta}-\theta_{0})^{\mathrm{T}}(\hat{\theta}-\theta_{0})<0.1$, otherwise we set $\theta_{0}=\hat{\theta}$ and re-iterate until the stopping criteria is satisfied.

B.3 Approximation of the Fisher Information

Louis (1982) and Oakes (1999) provide formulas for the Fisher information of the observed likelihood. We follow the recommendation of McLachlan and Krishnan (2007), arguing that Louis’s formula is best suited for the MCEM and provides a conservative measure of the standard errors. Therefore, we calculate the observed information based on

[TABLE]

Note that the second term of (16) depends not only on the first and second but also on the third and fourth conditional moment of the truncated multivariate normal and cannot be evaluated analytically therefore.

In order to approximate the observed information we are using the results of Robert and Casella (2004, p. 187) and Kang et al. (2013, Section 3.7) that allow for an approximation of the observed information based on the score and Hessian of the complete likelihood. Hence, we can use the results from Section B.1 for the following procedure.

We draw $w=1\,000$ times potential realizations $\tilde{Y}_{s,sim}$ from the truncated version of

[TABLE]

using the package TruncatedNormal (Botev, 2017). Those are stored for each draw $s$ in a vector $\tilde{Y}^{*}_{s,sim}=(\tilde{Y}_{o},\tilde{Y}_{s,sim})$.

Then we calculate the score and Hessian from equations (12) and (13) $w$ times, where we replace $\tilde{Y}$ by $\tilde{Y}^{*}_{s,sim}$ in each equation and index them by $s$, allowing to calculate the empirical version of (16) by approximating the expectations by means.

[TABLE]

This gives an estimator for the observed information. Standard errors are obtained by the square root of the diagonal elements of the inverted approximated matrix.

B.4 Data Transformation

In Figure 10 we show the distribution of the observed log-transformed response variable. The data is pooled over all years and standardized to have mean zero and variance one. The panel on the left side shows a kernel density estimate and the panel on the right gives a Q-Q plot.

B.5 Conditional Probabilities

Based on the fitted coefficients $\hat{\theta}_{t}=(\hat{\beta}_{t},\hat{\bm{\rho}}_{t},\hat{\sigma}_{t}^{2})$ and our model assumptions, we can represent the joint distribution of the latent utility network $Z_{t}$ via a multivariate normal

[TABLE]

where $\hat{\mu}_{t}=B(\hat{\bm{\rho}_{t}})X\hat{\beta_{t}}$ ands $\hat{\Sigma}_{t}=B(\hat{\bm{\rho}_{t}})(B(\hat{\bm{\rho}_{t}}))^{T}\hat{\sigma}^{2}_{t}$. Given that, define $Z_{t,-ij}$ as the $(N-1)$-dimensional vector, containing all entries of $Z_{t}$ except $Z_{t,ij}$. Additionally, for example in the case that $ij$ is the first entry of $Z_{t}$, rearrange $\hat{\Sigma}_{t}$ such that

[TABLE]

Then, the conditional distribution of $Z_{t,ij}$ is given by a univariate normal distribution

[TABLE]

We are interested in a possible state of the network, where the latent utility is allowed to be greater $c_{t}$. Therefore, we insert the expectation for the non-observed utility in $Z_{t,-ij}$ and denote this by $\tilde{Z}_{t,-ij}$. Consequently, we can calculate the probability of $Z_{t,ij}$ being greater than $c_{t}$ using

[TABLE]

The probability $\pi_{t,ij}$ can be interpreted as the probability that the latent utility of a transfer from country $i$ to country $j$ is higher than the threshold $c_{t}$ conditional on the covariates $X_{t}$ and the remaining network, where no transfer is restricted to be smaller $c_{t}$.

B.6 Simulation study - Endogenous effects and approximation of censored variables

In order to analyse the properties of our estimator, we use the following Data Generating Process (DGP1)

[TABLE]

Here, $q_{0.75}(Z)$ denotes the 75% quantile and we are censoring the network towards an observed density of $0.25$. Note, that DGP1 is not subject to under-reporting and all censored responses are in fact below the censoring threshold. The results of running DGP1 $100$ times and applying the estimation procedure are summarized in Figure 11. On the left panel, we show the true but censored values against the expected values from the last E-Step, together with contour curves and a non-parametric fit for the mean in solid black. It can be seen that the expected values approximate the latent variables very well. The right panel of Figure 11 shows boxplots for the difference between the true values of $\bm{\rho}$ and the estimated parameters. It indicates that as we are able to find unbiased estimates of the endogenous parameters despite the enormous amount of censoring.