Markov Chain Models of Refugee Migration Data

TL;DR

This paper explores the use of Markov chain models to analyze refugee migration data, demonstrating improved data fit and efficiency over existing agent-based models in the context of the Burundi refugee crisis.

Contribution

It introduces a Markov chain approach for refugee migration modeling, showing its advantages over agent-based models in accuracy and computational efficiency.

Findings

Markov chain models better match refugee migration data.

The approach is more computationally efficient.

Application to Burundi data demonstrates practical utility.

Abstract

The application of Markov chains to modelling refugee crises is explored, focusing on local migration of individuals at the level of cities and days. As an explicit example we apply the Markov chains migration model developed here to UNHCR data on the Burundi refugee crisis. We compare our method to a state-of-the-art `agent-based' model of Burundi refugee movements, and highlight that Markov chain approaches presented here can improve the match to data while simultaneously being more algorithmically efficient.

Tables

Table 1. Table 1: All time Average ARD ⟨ E ⟩ delimited-⟨⟩ 𝐸 \langle E\rangle , and Clean Average ARDs ⟨ E ⟩ 30 subscript delimited-⟨⟩ 𝐸 30 \langle E\rangle_{30} and ⟨ E ⟩ 100 subscript delimited-⟨⟩ 𝐸 100 \langle E\rangle_{100} for each model.

	Initial Graph	Camp-Adjusted	Graph-Adjusted	Time-Adjusted
$⟨E⟩$	0.41	0.22	0.32	0.37
${⟨E⟩}_{30}$	0.38	0.21	0.28	0.34
${⟨E⟩}_{100}$	0.37	0.21	0.24	0.28

Table 2. Table 2: All time Average ARD ⟨ E ⟩ delimited-⟨⟩ 𝐸 \langle E\rangle , and Clean Average ARDs ⟨ E ⟩ 30 subscript delimited-⟨⟩ 𝐸 30 \langle E\rangle_{30} and ⟨ E ⟩ 100 subscript delimited-⟨⟩ 𝐸 100 \langle E\rangle_{100} for each model.

	Initial Graph	Camp-Adjusted	Flee
$⟨E⟩$	0.41	0.22	0.29
${⟨E⟩}_{30}$	0.38	0.21	0.27
${⟨E⟩}_{100}$	0.37	0.21	0.28

Equations

$$f(i,j,k)=\min (f(i,j,k-1),f(i,k,k-1)+f(j,k,k-1)),$$

$$\sum _{1\le j\le n,\text{leavevmode}i\ne j}{I}_{ji}(t)=1-R_i.$$

$$I_{ji}(t)=\frac{\alpha \cdot \beta }{w_{ij}},$$

$$B(t)=\{N_1(t),\text{leavevmode}{N}_2(t),\cdots ,N_n(t)\}\text{leavevmode}.$$

$$A(t)=\left[\begin{array}{cccc}{P}_{11}(t)& P_{12}(t)& ...& P_{1n}(t)\\ P_{21}(t)& P_{22}(t)& ...& P_{2n}(t)\\ ⋮& ⋮& ...& ⋮\\ P_{n1}(t)& P_{n2}(t)& ...& P_{nn}(t)\end{array}\right]\text{leavevmode}.$$

$$[B(t+1)]_{i}=\big{[}A(t)\cdot B(t)\big{]}_{i}=P_{i1}(t)N_{1}(t)+P_{i2}(t)N_{2}(t)+\dots+P_{in}(t)N_{n}(t)\leavevmode\nobreak\ .$$

$$B(t^{\prime })=A(t^{\prime }-1)\cdot A(t^{\prime }-2)\cdots A(0)\cdot B(0)\text{leavevmode}.$$

$$g_{i}(d,t)\big{|}_{d\geq D}=\{0,0,\dots,0,1,0,\dots,0\}\leavevmode\nobreak\ .$$

$$g_{i}(d,t)\big{|}_{d<D}=\sum_{1\leq j\leq n}I_{ji}(t)g(j,d+w_{ij}).$$

$$R_{i}=\left\{\begin{array}[]{cc}0&v_{i}\leavevmode\nobreak\ {\rm is\leavevmode\nobreak\ a\leavevmode\nobreak\ conflict\leavevmode\nobreak\ site}\\ 0.7&v_{i}\leavevmode\nobreak\ {\rm is\leavevmode\nobreak\ a\leavevmode\nobreak\ neutral\leavevmode\nobreak\ site}\\ 0.999&v_{i}\leavevmode\nobreak\ {\rm is\leavevmode\nobreak\ a\leavevmode\nobreak\ camp\leavevmode\nobreak\ site}\end{array}\right.\leavevmode\nobreak\ .$$

$$T$$

$$\displaystyle=\frac{1}{N}\Big{(}MN(M)-\big{(}N(M-1)+N(M-2)+\dots+N(0)\big{)}\Big{)}.$$

$$E_i(t)=\frac{|N_i(t)-{\pi }_i(t)|}{U(t)}.$$

$$E(t)=\sum _{x_i\in X}{E}_i(t).$$

$$⟨E⟩=\frac{1}{396}\sum _{j=0}^{395}{E}_j=\frac{1}{396}(E(0)+E(1)+\cdots +E(395))\text{leavevmode}.$$

$$⟨E{⟩}_z\equiv \frac{1}{396-z}\sum _{j=z}^{395}{E}_j\text{leavevmode}.$$

Full text

1 Introduction

A record 65.6 million refugees were forcibly displaced in 2017 [1]. The economic disparities between different parts of the world, and spontaneous violent ethnic or political upheavals, have made migration one of the most significant issues of the modern age. Notably, the civil war in Syria triggered large numbers of refugees destined for Europe or neighboring Arab nations; the ethnic cleansing of Rohingya Muslims in Myanmar led to many Rohingyas fleeing to adjacent nations such as Bangladesh; and armed conflict in Somalia, coupled with widespread drought and famine, scattered refugees throughout other nations in the Horn of Africa.

Accurate models of refugee movements would in principle allow us to predict the number of refugees who will arrive at a particular area or city, the date at which they will arrive, and the distribution of these refugees across multiple regions, a few days or even weeks in advance. This would allow refugee and governmental organizations to determine where to best distribute aid resources to maximize impact and efficiency.

Migration has been studied since the 19th Century along the lines of general heuristics and simple “gravity” models [2, 3, 4], inspired by Newton’s theory of gravity. More recently, a variety of algorithmic approaches have been proposed, e.g. [9, 5, 7, 6, 8, 10, 11], many based on agent-based Model approaches [13] which bring significant improvements in their predictive capabilities. In this work we outline an alternative approach to simulating modern refugee crises by modeling the migration of refugees via stochastic matrices or Markov processes (see e.g. [22, 23]). The models developed here more accurately predicts the local movements of refugees in near-real time and improves on existing models in the literature. To this end we have developed for modern refugee crises on a regional level in near-real time, applied it to real world data and demonstrated that it runs efficiently from an algorithmic standpoint.

This type of Markov chain modelling has been previously successfully applied to a variety of topics, a small selection include viral epidemics [24], the internet [25] (notably, Google’s PageRank algorithm [26] which determines the relevance of webpages), financial systems [27], and evolutionary biology [28]. Markov chain models have also been applied to study large-scale problems of long-distance and more global migration patterns [14, 15, 16, 17, 18, 19, 20, 21], but to our knowledge this presents the first application to local movements of people and in particular to the study of refugees. We shall argue that Markov chains are an efficient and powerful manner to model refugee movements and present an alternative to popular agent-based models.

This work is structured as follows: In Section 2, we outline the geographical component of our approach to modelling refugee movements. In Section 3, we outline the algorithmic component of the Markov chain migration model. In Sections 4 & 5, we apply several versions of the MC Migration Model to the recent refugee crisis in Burundi and present results. In Section 6, we compare our simulations with results from an existing model of the Burundi refugee crisis. In Section 7 we provide a summary and some concluding remarks.

2 Graphs from Geography

The first key component in modelling refugee movements is accurately encapsulating geographical information from the region of interest. The most common approach is to construct a weighted graph $G$, where the vertices $\{v_{1},v_{2},\dots,v_{n}\}\in G$ represent cities in the region and the edges in the graph represent roads between cities, weighted by the lengths of the roads. However, this leaves some degree of freedom in choosing which roads and cities to include, and how to assign weights; moreover refugee crises often occur in underdeveloped nations where population centers are small and the geographic data incomplete. For instance, one approach to determining which urban centers are significant enough to be included as vertices in the graph when reliable population data is unavailable is to utilizes nighttime images from satellites along with population estimates [6].

A geographic map is a projection of some physical region to $\mathbb{R}^{2}$ and cities can be represented as points in $\{v_{1},v_{2},\dots,v_{n}\}\in\mathbb{R}^{2}$. Subsequently, we construct a planar graph $G_{P}$ by identifying the vertices as $\{v_{1},v_{2},\dots,v_{n}\}$ and introducing edges which correspond to the most major roads in the region in such a manner that any two edges intersect at their endpoints only. As a result, edges typically only link cities with other cities in their immediate vicinity.

An edge connecting two cities $v_{i}$ and $v_{j}$ is assigned a weight $w_{ij}$ corresponding to the physical distance in kilometers between $v_{i}$ and $v_{j}$, as determined via Google Maps. The graph construction thus far follows that of a “Local Interaction Model” detailed in [6], however, since many refugees will move along minor roads not included in this most basic planar graph, we will refine this initial graph construction with certain heuristics detailed below.

For any two cities $v_{i},v_{j}$, let $d(i,j)$ be the minimal distance between $v_{i},v_{j}$ using edges from $G_{P}$, meaning $d(i,j)$ is the length of the shortest path from $v_{i}$ to $v_{j}$. We propose that there is some characteristic maximum distance $D$ that a refugee can travel in one day. The value of $D$ must be assigned with intuitive reasoning, rather than some prescribed method, and should incorporate information such as the geography of the region and the quality of roads. The value of $D$ is significant because it allows one to distinguish between different refugee movements based on their distance. We construct a refined graph $G$ of a given geographical region with cities $\{v_{1},v_{2},\dots,v_{n}\}$ by adding additional edges to $G_{P}$. Specifically, for each non-adjacent pair $v_{i},v_{j}\in G_{P}$ with $D\leq d(i,j)\leq 2D$, an additional edge $\{v_{i},v_{j}\}$ with weight $w_{ij}=D$ is introduced in the planar graph $G_{P}$. The resulting construction $G$ is typically non-planar.

The Floyd-Warshall algorithm [30, 31, 32] can be utilized to efficiently construct this new graph. The Floyd-Warshall algorithm takes a weighted graph with vertices $v_{1},v_{2},\dots,v_{n}$ such that each pair of vertices $v_{i},v_{j}$ is connected by an edge $e_{ij}$ of weight $w_{ij}$; if $v_{i},v_{j}$ are not adjacent, we adopt the convention $w_{ij}=\infty$. The algorithm then computes the shortest paths between every pair of vertices $v_{i},v_{j}$. It accomplishes this by letting $f(i,j,k)$ be the minimum length of a path between vertices $v_{i},v_{j}$ such that any intermediate vertices $v_{l}$ satisfy $l\leq k$. Then $f(i,j,0)$ is $w_{ij}$ because there can be no intermediate vertices and one can recursively calculate

[TABLE]

because the first term is the shortest valid path from $v_{i}$ to $v_{j}$ which does not pass through $v_{k}$, while the second term is the shortest valid path from from $v_{i}$ to $v_{k}$ and then from $v_{k}$ to $v_{j}$. The recursive step takes $O(1)$ time to compute each new value $f(i,j,k)$ for $1\leq i,j,k\leq n$ in terms of previously stored values $f(i,k,k-1),f(j,k,k-1)$, so the entire algorithm runs in $O(n^{3})$ time. Then $d(i,j)$ is simply $f(i,j,n)$, so all minimal distances are computed. To our knowledge, this approach of distinguishing between different pairs of vertices using their minimal distance in a planar graph has not been applied to migration before.

3 Markov Processes and Refugee Movements

In what follows we will combine the graph constructions outlined in the previous section with novel modelling approaches which we discuss next in order to generate a new class of migration models, whose application we explore in Section 4. Consider a graph $G$ with vertices $v_{1},v_{2},\dots,v_{n}$ such that some pairs of vertices $v_{i},v_{j}$ are connected by edges $e_{ij}$ of weight $w_{ij}$. We demonstrated our approach to constructing such a graph which represents a geographical region in Section 2 and we next explain how migration can be simulated on such a graph. To do this we adapt some ideas and conventions developed in previous approaches to modelling refugee movements using agent-based models, due to Groen [5].

The vertex in the graph indicate specific geographical locations. To model refugee migration it is useful to differentiate these nodes between three types of location, specifically

• •

Camps: Refugees in a refugee camp have a high probability of remaining in the camp, and refugees near a refugee camp have a high probability of migrating into that camp.

• •

Conflict: Refugees in a conflict site have a high probability of leaving, and refugees near a conflict site have a low probability of entering the conflict site.

• •

Neutral: A neutral vertex is a location which is neither a refugee camp or conflict site. Refugees in a neutral city have a moderate probability of entering or leaving.

A refugee camp is a location set up to house refugees and a conflict site is a city which has been the location of an armed conflict or other disaster. In principle a location can change type during the course of the simulation, for instance, as new refugee camps are established. Each vertex in our model carries this ‘type’ information as a label, along with three other labels:

• i).

Location name;

• ii).

location type (Camp/Conflict/Neutral);

• iii).

Time stamp $t\in\mathbb{N}$;

• iv).

Refugee population at a given time stamp $N(t)\in\mathbb{N}$.

To simulate the movement of refugees the population value $N(t)$ of each vertex evolves between each time stamp. Our model iterates through a series of discrete timestamps $t=0,1,2,\dots$ with each timestamp representing a single day and $t=0$ corresponding to the initial configuration of the system. The premise of our model is that on each day, refugees make a series of moves between locations in a region; once the total distance a refugee has travelled in a single step reaches or exceeds the maximum threshold $D$, the refugee cannot make any more movements on that day. Between timestamps the population of each vertex changes due to the movements of the refugees.

The path that a refugee takes through the system is determined by the maximum threshold $D$ and also the properties of the vertices. Specifically, in our model a refugee at vertex $v_{i}$ has a probability of $R_{i}$ of remaining in that state. We assume that $R_{i}$ depends only on status of $v_{i}$ as either a camp, conflict, or neutral site. In addition to this we next define the probabilities that a refugee at a given location migrates to a different location between timestamps. Consider two consecutive timestamps $t^{\prime}$ and $t^{\prime}+1$ and let $v_{i}$ and $v_{j}$ represent two adjacent vertices in $G$. For a refugee at vertex $v_{i}$ at timestamp $t^{\prime}$ which is able to migrate we denote by $I_{ji}(t^{\prime})$ the probability that a refugee at $v_{i}$ for $t=t^{\prime}$ moves to $v_{j}$ at $t=t^{\prime}+1$. We refer to $I_{ji}(t)$ as the intermediate probability at timestamp $t$.

Note that $I_{ii}(t)$ represents the probability that the refugee remains in city $v_{i}$, so $I_{ii}(t)=R_{i}$. Furthermore, we assume that a negligible number of refugees leave the system of urban centers, so the laws of probability imply that $I_{1i}(t)+I_{2i}(t)+\dots+I_{ni}(t)=1$ for each $1\leq i\leq n$ where $n$ is the number of vertices in $G$. It follows that

[TABLE]

In our model the probability that a refugee will move to a particular location is a function of the edge weight $w_{ij}$ (which encodes the geographical distances) and a scaling factor $\alpha$ that takes different constant values depending on whether the location is a camp, conflict or neutral site (we specify the values of $\alpha$ in Section 4 for a specific refugee crisis). For each vertex $v_{j}$ adjacent to $v_{i}$ we define the intermediate probability to be

[TABLE]

where $\beta$ is a constant scaling factor fixed such that the probabilities add to unity after the type rescaling, thus eq. (2) is satisfied.

This setup has the following desired properties:

• •

The total sum of the intermediate probabilities from any city $v_{i}$ is $1$.

• •

Refugees are more likely to migrate to refugee camps than to neutral cities.

• •

Refugees are more likely to migrate to neutral cities than to conflict sites.

• •

Refugees cannot directly migrate between non-adjacent cities.

• •

The probability of moving to an adjacent city $v_{j}$ is inversely proportional to the distance $w_{ij}$ separating the cities.

Note that this last point is in accordance with Ravenstein’s Laws of Migration [2], which state that migrants prefer to move shorter distances.

Thus far we have discussed how to encode the local preferences for refugee movements which is characterized by the intermediate probability defined above. We next introduce the notion of a transition probability $P_{ji}(t)$, being the probability that a refugee in city $v_{i}$ at the start of timestamp $t$, will move to city $v_{j}$ at the start of the next timestamp $t+1$.

Our approach maintains that a refugee may make multiple moves between cities in a single day until the refugee’s total distance travelled exceeds $D$. Each individual move is made according to the intermediate probabilities $I_{ji}$, while the transition probabilities $P_{ji}$ measure the net effect of these movements taken in series. Thus the transition probabilities are completely determined by the choice of intermediate probabilities.

Suppose that at timestamp $t$ that the refugee population of city $v_{i}$ is $N_{i}(t)$ for $1\leq i\leq n$, we define the $n$ dimensional ‘population’ vector which is a function of the timestamp $t\in\mathbb{N}$

[TABLE]

For a system with transition probabilities $P_{ij}(t)$ for $1\leq i,j\leq n$, the stochastic matrix of the system at time $t$ can be denoted as

[TABLE]

The entries of the matrix are each probabilities with the rows summing to one, thus the matrix describes transitions in a Markov chain.

The stochastic matrix $A(t)$ dictates the probability of moving from $v_{j}$ to $v_{i}$ between timestamp $t$ and $t+1$. Thus the expected population vector at timestamp $t+1$ can be described compactly as $B(t+1)=A(t)B(t)$. Moreover, the $i$th entry of $A(t)B(t)$ gives the expected number of refugees at vertex $v_{i}$ at timestamp $t+1$ given by

[TABLE]

More generally for a given timestamp $t^{\prime}$ the vector $B(t^{\prime})$ can be defined from the initial population vector $B(0)$ recursively

[TABLE]

Thus at the start of each sequential timestamp the vertex populations are inherited from the final populations at the end of all operations of the previous timestamp. This is traditional Markov process of an iteratively evolving system. More specifically, since the form of $A$ is not fixed but intermittently evolves over the timestamps if a neutral city becomes a camp or conflict site, this is a non-homogeneous Markov chain [23] .

This Markov system outlined above is advantageous for modelling migration because its allows for a relatively compact model solely in terms of the transition probabilities $P_{ji}(t)$ and concerns about intermediate probabilities and about refugees moving across multiple cities in one day are irrelevant because they are already encapsulated into the transition probabilities. As a result, the main difficulty in implementing a stochastic matrix model is to compute the values of $P_{ji}(t)$. Our algorithm for evolving the vertex populations between two timestamps performs a set of order operations and thus one can break the problem of computing transition probabilities into smaller pieces. This permits the application of dynamic programming techniques [29] which are essentially an efficient form of recursion. To implement this we introduce a new location specific $n$-tuple $g_{i}(d,t)=\{s_{1}(d,t),s_{2}(d,t),\dots,s_{n}(d,t)\}_{i}$ associated to vertex $v_{i}$. The elements $s_{j}(d,t)$ of $g_{i}$ represent the probability that the refugee will be at vertex $v_{j}$ at the end of timestamp $t$, which depends on the distance $d$ a refugee has previously travelled during the current timestamp $t$ such that they currently reside in city $v_{i}$.

Specifically, if $d\geq D$, then the refugee cannot make any further movements to different cities, so will remain in city $v_{i}$ for the remainder of the timestamp. Therefore

[TABLE]

where the single $1$ occurs in the $i$th entry of the vector. If $d<D$, then for each $1\leq j\leq n$, the refugee moves from $v_{i}$ to $v_{j}$ with intermediate probability $I_{ji}(t)$; at the end of this move, the refugee is now in city $v_{j}$ and has moved a total distance of $d+w_{ij}$. At this point, the probability distribution for the final destination of the refugee is $g(j,d+w_{ij})$. It follows that

[TABLE]

Note that when $v_{i},v_{j}$ are not adjacent, our sum contains a term $I_{ji}(t)g(j,d+w_{ij})=0\cdot g(j,\infty)$; we adopt the convention $g(j,\infty)=\{0,0,\dots,0\}$ for all $j$ to ensure the expressions are well-defined, though the multiplication by zero implies the exact value of $g(j,\infty)$ is irrelevant.

Thus eq. (9) yields a recursive formula for $g_{i}(d,t)$ in terms of expressions of the form $g_{j}(d^{\prime},t)$, where $d^{\prime}=d+w_{ij}>d$. As a result, we can efficiently compute all vectors $g_{i}(d,t)$ by first computing all vectors of the form $g_{i}(D-1,t)$ and storing the results, then computing all vectors of the form $g_{i}(D-2,t)$, and so on; eq. (8) handles terms $g_{i}(d,t)$ with $d\geq D$ that occur in the recursion. Computing each $g_{i}(d,t)$ in terms of $n$ terms of the form $I_{ji}(t)g_{j}(t)(d+w_{ij},t)$ takes $O(n^{2})$ time: It is essentially summing $n$ vectors of the form $g_{j}(d+w_{ij},t)$, each of which has already been computed and is $n$ dimensional. Therefore, this algorithm computes the $nD$ vectors of the form $g_{i}(d,t)$ for $1\leq i\leq n,d<D$ in $nD\cdot O(n^{2})=O(n^{3}D)$ time.

The transition probability $P_{ji}(t)$ is equivalent to the probability that a refugee currently located at city $v_{i}$ during timestamp $t$ who has travelled $d=0$ distance ends up at city $v_{j}$ by the end of timestamp $t$. As a result, $P_{ji}(t)$ is the $j$th entry in $g_{i}(0,t)$ for all $1\leq i,j\leq n$, so the computation of all $g_{i}(d,t)$ terms allows one to compute all the transition probabilities $P_{ji}$(t).

The transition probabilities are constant over time unless the intermediate probabilities $I_{ji}(t)$ change which only happens when a city becomes a refugee camp or conflict site and whilst this occurs in our models, it is not a common occurrence. As a result, the above procedure should only run a few times during each refugee crisis. The computation of transition probabilities from intermediate probabilities completes the algorithmic aspect of our model.

4 Application: Modelling the Burundi Crisis

We utilize data on the Burundi crisis from a previous model of the crisis called Flee due to Suleimenova-Bell-Groen (SBG) [8], who in turn obtained their data from the United Nations High Commissioner for Refugees (UNHCR). The UNHCR data extracted by SBG [8] runs from May of 2015 to June of 2016, for a total of $396$ days, so the timestamps in the model range from $0\leq t\leq 395$. There were five active refugee camps in the Burundi crisis, each established by the UN at the cities of Nyarugusu, Nduta, Nakivale, Mahama, and Lusenda. Figure 1a shows the 30 major cities in and around Burundi which were modelled in Flee and Figure 1b illustrates the graph used to model the system constructed and used in SBG [8]. Figure 1c shows a non-planar graph representing additional routes connecting relevant locations, which we constructed by adding further edges to the graph $G_{P}$ in panel (b), as detailed in Section 3.

In SBG [8] the authors assumed that the probabilities $R_{i}$ of remaining in a city $v_{i}$ are

[TABLE]

These estimates seem reasonable so we maintain them in our model too. We also take a list of major armed conflicts in Burundi from SBG’s Flee model [8] and turn the corresponding locations into conflict sites at the appropriate timestamp. The first such conflict took place on May 1st, 2015, which we set as the beginning date $t=0$ of the model. SBG [8] estimate that refugees travel a maximum distance of $200$ km in one day.222The supplementary materials of [8] presents sensitivity testing for Flee and shows that the results are robust under modest variations of the maximum distance $D$ and the staying probabilities $R_{i}$. Arguably, this overly optimistic as most refugees travel on foot [8], so we instead set a daily distance threshold of $D=120$ (in units of km), corresponding to $15\text{ hours}$ at 8 km/h.

Next, we initialize the population vectors $B(t)$ to match the number of refugees in the Burundi region; again we obtained this data from [8], who retrieve it from UNHCR refugee camp registration data. SBG [8] notes that it is incredibly difficult to predict the number of refugees created by armed conflicts, but rough estimates can be obtained through various means such as satellite images [6]. As a result, rather than predicting the number of refugees resulting from a crisis, SBG [8] aim to predict the distribution of these refugees given the total number of them; our model also focuses on this objective. We will return in Section 5 to directly compare our results with those from the Flee model of SBG [8].

Since $B(t+1)=A(t)B(t)$ represents the expected distribution at time $t+1$ of the refugees from time $t$, the total number of refugees in $B(t+1)$ is the same as the total number of refugees in $B(t)$. This implies that any increases in the number of refugees are not reflected in the matrix multiplication and therefore must be performed through an additional process. Specifically, we denote by $U(t)$ the number of refugees in Burundi refugee camps at time $t$ according to the UNHCR. After each timestamp $t$, the population discrepancy $U(t+1)-U(t)$, between timestamp $t$ and timestamp $t+1$, is compute and then an additional $U(t+1)-U(t)$ refugees is randomly added into the population vector $B(t)$ by distributing them into the conflict sites in the model with probability proportional to the population of the conflict site. This ensures the model always contains an accurate number of refugees.

Due to this mechanism of adding refugees, after timestamp $t$, the total number of refugees in the model is the same as the total number of refugees in UNHCR camps after day $t$. As a result the model continually under populates the graph and therefore under predicts the number of refugees in each camp, because many refugees in the model are still making their way to camps. To account for this delay between the creation of refugees and their arrival in camps in the model, when the refugee populations in each refugee camp are extracted, following [8] we rescale each such population upwards by a constant factor so that the total number of refugees in the camps matches the UNHCR data for each day. This rescaling occurs in three of the four Markov chain models we outline below and in Flee [8]. Because our stated goal is to predict the distribution of refugees in camps given the total number of refugees, this rescaling is unimportant. However, since this scaling is somewhat ad hoc one of the Markov chain variants we outline below is designed to remove the need for this rescaling.

We now apply the modelling algorithms outlined in Sections 2 & 3 to the recent Burundi refugee crisis, due to a civil war, focusing on a set data on camp registrations during a period from 2015-2016 . We run four different models using Markov chains, each of which is a variation on the principles outlined in Sections 2 & 3, as we outline below:

Markov Chain: Initial graph. We apply the Markov Chain approach to the planar graph $G_{P}$ constructed in [8] from an analysis of major roads in Burundi on Bing Maps, see Figure 1b.

Markov Chain: Graph-adjusted. We additionally apply our Markov chain approach to the modified graph $G$ which we derive from $G_{P}$ using the procedure from Section 2, show in Figure 1c, in what we call the graph-adjusted model.

Markov Chain: Camp-adjusted. Long movements to refugee camps are systemically assigned smaller intermediate probabilities in our model. This heuristic makes sense for movements to neutral cities because refugees prefer making a series of shorter movements rather than a single long movement [2], but this reasoning breaks for refugee camps because they are refugees’ final destinations rather than intermediate locations. In addition, there are a variety of motivations for refugees to prefer faraway refugee camps, such as perceived safety and distance from danger. To realise this observation, we create an additional graph $G^{\prime}$ by taking edges $e_{ij}$ in $G$ that represent long movements to refugee camps and reassigning these edges with weight $w_{ig}=D$ so that long movements to refugee camps have the same chance of occurring as medium movements to camps. We call this variation the camp-adjusted model.

Markov Chain: Time-adjusted. We now present a fourth Markov chain model, which we call the time-adjusted model, which removes the need for the overall rescaling to match the total refugee population at the end of each timestamp. Denote by $T$ the average time it takes for a refugee in the model to travel from a conflict site to a refugee camp. Then if there are $U(t+T)$ refugees in the model at timestamp $t$, most of these refugees should arrive in refugee camps within $T$ days, so the number of refugees residing in camps within the model will be approximately $U(t+T)$ at timestamp $t+T$. We can use this observation to remove the need for the rescaling at the end of each timestamp by introducing an appropriate number of refugees in the model $T$ days before they are observed. Specifically, at timestamp $t=0$ we introduce $U(T)$ refugees instead of $U(1)$ refugees, and at the end of each timestamp $t$ we add an additional $U(t+T+1)-U(t+T)$ refugees rather than $U(t+1)-U(t)$ refugees.

It remains to compute $T$, the average travel time for a refugee from a conflict site to a refugee camp, which can be found by running a mock simulation where the entire population $N$ of the conflict sites become refugees at timestamp $t=0$, and then let $N(\delta)$ denote the total camp population after $\delta$ days for $1\leq\delta\leq M$ for some suitably large constant $M$. Then $N(\delta+1)-N(i)$ refugees take exactly $\delta+1$ days to arrive, so the average travel time is

[TABLE]

For $M=100$ we find $T=16$ and thus this is the value we use in the time-adjusted model.

5 Results: Modelling the Burundi Crisis

We will now present the results of the Burundi models outlined above. Figure 2 shows five different graphs which are model results for number of refugees in the major regional camp as a function of time. Each graph represents one of the five different camps relevant to the Burundi crisis, and in each graph we plot the refugee population in the given camp as a function of the timestamp (days). Within each graph we show five curves: the four solid curves correspond to the four Markov chain models outlined in Section 4, the dashed curve is the UNHCR camp population data, and the shaded region indicates a 10% error on the data.

It is difficult to determine from these population plots alone how effective each model is at predicting refugee populations in the five camps. To quantify the goodness of fit we now introduce metrics to measure the success of a given refugee model. Denote by $X$ the set $\{\text{Nyarugusu},\text{Nduta},\text{Nakivale},\text{Mahama},\text{Lusenda}\}$ of relevant refugee camps. For each camp $x_{i}\in X$, let $N_{i}(t)$ be the population of that camp at timestamp $t$ in a given model, and let $\pi_{i}(t)$ be the population of that camp at day $t$ according to UNHCR data. Note that $U(t)$ is precisely the sum of $\pi_{i}$ over each of the camps: $\sum_{x_{i}\in X}\pi_{i}(t)=U(t)$. The first measure we consider is the Averaged Relative Difference (ARD) for a given camp $x_{i}$ on day $t$, given by

[TABLE]

The total ARD on a given timestamp $t$ is then the sum over the ARD values of all of the camps

[TABLE]

Moreover, it is also useful to consider time averaged ARDs for a given model, by averaging over all the daily ARD values in a certain period. In particular, we consider the Average ARD $\langle E\rangle$ over the entire period covered by the data (up to $t=396$)

[TABLE]

Figure 3 shows the total ARD $E(t)$ for each model as a function of time. This provides a good comparison tool to assess how well the models match the data. Observe that the Markov chain model using the ‘Initial graph’ (Figure 1b) consistently has relatively high error margins, while the graph-adjusted and camp-adjusted models (Figure 1c) are substantially improved.

Notably, there is potentially a lot of noise in the first few days of every model and in actual refugee crises, thus we expect the quality of the data to improve over time, as registration practices improve and the situation in the camps becomes more stable. To better analyze the long-term accuracy of the models we also consider the averaged ARD omitting the earliest period of the data, since this likely provides a more reliable measure of ‘goodness of fit’. Specifically, we call the Average ARD omitted the first $z$ days (counting from day $t=0$) the $z$-day Clean Averaged ARD and denote this by $\langle E\rangle_{z}$ with

[TABLE]

We will consider the 30-day and 100-day Clean Averaged ARD, denoted $\langle E\rangle_{30}$ and $\langle E\rangle_{100}$, the values of which are shown in Table 1 for the four Markov chain models of Section 4.

From the statistics from Table 1 we can draw some general conclusions about our models. Firstly, the Camp-Adjusted Model has a lower Average ARD than every other model in each of the three numerical measures $\langle E\rangle$, $\langle E\rangle_{30}$, and $\langle E\rangle_{100}$. Furthermore, the results clearly demonstrate the merit of the geographical adjustments introduced in Section 2, with the Camp-Adjusted and Graph-Adjusted Models yielding the lowest long-term average ARDs. That the Camp-Adjusted and Graph-Adjusted models give an improved match to data compared to graphs which represent the major road system may be indicative of refugees travelling off-road, or on minor roads not well documented on Bing or Google Maps. It would be interesting to further validate this conclusion over additional conflicts and longer term data sets.

The Time-Adjusted Model begins with relatively inaccurate results, but the average error quickly drops, as expected because this is a time-sensitive model. However, omitting the first 100 days and considering $\langle E\rangle_{100}$, the Time-Adjusted Model gives a comparably low average ARD to the other models, indicating that we were successfully in removing the rescaling of refugee camp populations at the end of each timestamp while still producing a relatively accurate model. This is significant because removing the rescaling process results in a more transparent model and less reliance on UNHCR data.

A key parameter in the Markov chain models is the distance threshold $D$ which determines the maximum distance a refugee can travel in one day. This parameter is also a source of some uncertainty since, whilst our assumed choice $D=120$ (in units of km) seems quite reasonable, it could conceivably vary by moderate margins depending on road quality, climate, and transportation options. Moreover, for the Camp Adjusted model as $D$ is varied the number of edges in the graph varies as described in Section 2. Thus it is prudent to explore how varying $D$ impacts the late time distribution of refugees. We show the impact of varying $D$ on the camp populations at $t=396$ (the final entry in the data set) in Figure 4, the solid curves give the prediction for the Markov Chain: Camp Adjusted model, and the dashed lines indicate the value reported in the UNHCR data. Observe that for $D\ll 100$ then the fit to the late time data becomes marginally worse, however, for $D>100$ the results are relatively insensitive to changes in $D$ which provides some confidence that we need not be concerned by the exact value or indeed the possibility local variations in the maximum travel distance.

6 Discussion: Comparison to Agent-Based Models

We next discuss the differences between our Markov chain model and previous approaches using agent-based models, with specific reference to one state-of-the-art example due to Suleimenova, Bell, & Groen [8] called Flee. We will briefly summarize the main components of Flee and then compare our results to those obtained from running Flee (which is publicly available).

Flee [8] implements a drastically different approach known as agent-based modelling [13]. An agent-based model consists of a series of fixed states and a group of agents. Individual agents transition through these states according to the rules of the system, updating their current state at each integer timestamp. In modelling migration, a common approach is to define a weighted graph to contain the geographic information of a region of interest, and introducing one agent to represent each individual the region. These agents then migrate between the vertices of the graph according to some predefined algorithm. This approach has been well-studied and utilised in areas such as disaster-driven migration [34], the ongoing conflict in Syria [35], and in [8] applied to conflicts in Burundi, Mali, and the Central African Republic.

As in our Markov chain model, Flee [8] consists of a graph storing geographic information and a prespecified algorithm used by refugees to move between cities of the graph. We noted earlier that we derive the planar graph (Figure 1b) in our Markov chain model from Flee. Meanwhile, their agent-based model differs significantly from our Markov chain model, because Flee treats each refugee in the simulation as an entirely separate entity making independent choices. Whereas we model refugee movements in terms of general probabilities of transitioning between states, Flee iterates over all agents in the system during every timestamp and simulates for each individual agent its expected movement through various cities.

We noted earlier that the heuristic we use for properly initializing the population vectors at any given timestamp is derived from Flee’s heuristic for generating new agents at each day of the simulation, the number of new agents added to the system is computed using real data about refugee numbers. This is acceptable because Flee aims to predict the end distribution of refugees over a network rather than the actual number of refugees displaced, as such, a rescaling of the distribution is performed at the end of every timestamp to ensure the total number of refugees in the camps in Flee matches the corresponding total according to UNHCR data, the same rescaling we use in three of the Markov chain models.

Despite some similarities between Flee and our Markov chain models, differences in algorithms lead to significantly different results. We now compare the evolution of the camp populations and the ARD measures defined in eqns. (11)-(13) to quantitively compare our Markov chain models to Flee. Figure 5 is analogous to Figure 2, but illustrates instead the change in the populations of each camp as predicted by Flee (which match results presented in [8]) to our Markov chain models implemented on the Initial Graph (using same graph as [8]) and the camp adjusted graph (Figure 1c).

Additionally, in Figure 6 we show the total ARD as a function of timestamp for the three models displayed in Figure 5, which gives provides a quantitive assessment of the relative goodness of fit of each model. It is noteworthy that by the end of the simulation, each of the five models has more or less converged to a constant value, and Flee has the highest total ARD. The trends are captured concisely in Table 2 which shows the all time Average ARD $\langle E\rangle$, and the Clean Average ARDs $\langle E\rangle_{30}$ and $\langle E\rangle_{100}$ for each model.

The Camp-Adjusted Model and Flee display relatively stable average ARDs of approximately $0.21$ and $0.28$ over time, while the other models take longer to converge to lower average ARD values. In particular, note that in terms of long-term predictive power, which is the main focus of refugee models, the Markov chain: Camp Adjusted model outperforms both Flee and the Markov Chain: Initial Graph. The Camp-Adjusted model has a value of $\langle E\rangle_{100}$ for the is roughly $76\%$ that of Flee, implying a $24\%$ reduction in the long-term error.

In [8] to quantify their results the authors use instead the camp-wise ARD $E_{i}(t)$ for a given day as the main metric of success in modelling the Burundi crisis. Here we have focused instead on the total ARD which we believe is a clearer indicator of overall modelling accuracy, however the camp-wise ARDs for each model are provided in the Appendix to provide a more direct comparison and one can draw similar conclusions on the goodness of fit of each model from looking at the ensemble of the individual camp ARDs.

Thus we conclude that Markov chain models can give a modest improvement on the goodness of fit compared to Flee, moreover, Markov chains yield a significantly more efficient simulation tool than an agent-based approach in terms of runtime and code length. In terms of code length, the components necessary to make complex agent-based models such as Flee operate are spread over multiple folders, dozens of files, and thousands of lines in code [8]. In contradistinction, each of our Markov chain models is contained within a single file and less than 500 lines of code. While program length is no measure of efficiency or suitability, the maintenance costs of code and the likelihood of bugs are directly correlated with program length.

To quantify the complexity of the agent-based Flee, consider a total of $R$ refugees and $n$ cities in a region of interest, for a model which runs $t$ time steps. At each timestamp, Flee iterates through all agents in the system to update their status, because agent-based models treat each agent separately. Furthermore, each agent runs through the list of $n$ cities before deciding which city to migrate to, thus each of $R$ agents takes $O(n)$ time to update at each timestamp, leading to a total time complexity of at least $O(nRt)$.

Analyzing the overall time complexity of our Markov chain models is slightly harder. If there are $n$ cities in the region of interest, then we noted in Section 2 that our graph modification algorithms utilizing Floyd-Warshall run in $O(n^{3})$ time, however, these algorithms only run once, at the very beginning of the simulation, to establish the graph. Meanwhile, we noted in Section 3 that we compute the transition probabilities $P_{ji}$ in $O(n^{3}D)$ time. Once the transition probabilities are computed once, they do not need to be recomputed, as we reuse the same stochastic matrix $A(t)$ over multiple timestamps, unless a fundamental change is made to the cities in the graph, for instance, a neutral city turning into a conflict site changes the intermediate probabilities and therefore the transition probabilities of the system. The number of such updates to vertices is on the order of the number of vertices $n$, because neutral cities can become conflict sites, but the opposite is much rarer. As a result, computing transition probabilities takes around $O(n)\cdot O(n^{3}D)=O(n^{4}D)$ time. Finally, at each timestamp, the Markov chain only performs a simple multiplication of an $n\times n$ matrix with an $n$-dimensional vector, which only takes $O(n^{2})$ operations, thus updates to the system comprise a total of $O(n^{2}t)$ time. Thus we estimate the time complexity of the Markov chain models is $O(n^{3}+n^{4}D+n^{2}t)$.

In the case of Burundi, $R\sim 10^{5}$ refugees, $t=396\sim 10^{2}$ days, and there are $n\sim 10$ vertices, while $D\sim 10^{2}$, thus Flee runs in $O(nRt)=O(10^{8})$ operations, while the Markov chain models are $O(n^{3}+n^{4}D+n^{2}t)=O(10^{6})$ operations, because the $n^{4}D$ term dominates the asymptotic. Indeed, this order of magnitude difference in runtimes is apparent when executing these different simulations on a standard desktop or laptop. In general, migration models will typically only involve significant urban centers, so the value of $n$ is reasonably small, while $R$ can be incredibly large, making agent-based model less feasible than Markov chains implementations.

7 Concluding Remarks

In this work we have implemented Markov chain models of refugee migration in refugee crises based on new heuristics about the connections between refugee movements and local geography. By applying our model to the Burundi refugee crisis and comparing with an existing agent-based model of that same crisis, we concluded that our approach using Markov chains was more efficient, reduced unnecessary complications in agent-based modelling, and exhibited a $24\%$ reduction in long-term prediction errors. The agent-based model Flee was explicitly validated against three separate conflicts (Burundi, Mali, and the Central African Republic) in [8], whereas for brevity we have restricted our analysis here to just the Burundi conflict. However, we highlight in passing that the Markov Chain models developed here similarly provide good fits (comparable to Flee) to these other two conflicts when applied to the data and graphs of Mali and the Central African Republic presented in [8].

It should be noted that collecting data on refugee registrations at camps is a challenging task and it is highly feasible that there could be sizeable discrepancies between the UNHCR data and real camp populations. Moreover, refugees that leave a camp for a new destination, say another camp or perhaps a relatives house in a neutral city, are likely not to see deregistering as a priority and thus the data my not accurately reflect departures from a camp. We have endeavored to account for these potential discrepancies in the data by including a 10% error margin in Figures 2 & 5, however this is a rather blunt heuristic and a better understanding of the intrinsic errors in the data would be highly desirable. Applying our model to higher quality and larger data sets would allow us to better refine the Markov chain models presented here.

As a result of the recent increases in refugee crises, there are many reasons as to why improved refugee models would be useful. Accurate refugee models would allow us to predict the number of refugees who will arrive at a particular area or city, the date at which they will arrive, and the distribution of these refugees across multiple regions, a few days or even weeks in advance, given good information about the causes of the migration. This would be extremely beneficial as accurate modelling could allow refugee and governmental organizations to determine where to best distribute aid resources to maximize impact and efficiency. Predictive models would furthermore allow cities and regions to take the appropriate measures to accommodate a large influx of people seeking shelter, which can typically lead to highly disruptive situations. Indeed, with the onset of global climate change, increasing inequality, and decreasing global stability it is highly likely that the number of displaced peoples will increase significantly in the future [36]. Thus there is a growing need for better simulations and in this work we have advocated for adopting Markov chain based models in order to more accurately and efficiently model these crises.

Acknowledgements

We thank William Jones, Laura Schaposnik, and especially Derek Groen for invaluable discussions, as well as Tanya Khovanova and Claude Eicher for comments on the draft. JU is grateful to the Simons Center for Geometry and Physics (Program: Geometry and Physics of Hitchin Systems) and New College, Oxford for their hospitality and support. This work was completed as part of the MIT PRIMES Program.

Appendix: ARDs for Individual Refugee Camps

In this appendix we present the ARD for each camp as a function of timestamp. The ARD for a camp gives a good measure for the amount of error present in modelling a single camp as a fraction of total refugee population. Note that Figures 3 & 6, which show the total ARD for each model on each day essentially “sum” of the five graphs in Figure 7.

Figure 7 shows the ARDs in each of the five refugee camps of each of the five models, the four Markov chain models developed here, and SBG’s Flee model [8]. Observe that in two of the camps, Lusenda and Nakivale, all five models perform quite well, leading to an almost indistinguishable long-term ARD of approximately zero. However, this trend does not hold in the other three camps.

In Nduta, all three models utilizing our graph modification algorithm from Section 2 perform significantly better than than Flee and the Markov Chain Model: Initial graph, which both use the unmodified graph. In Nyarugusu, Flee initially produces relatively accurate results, but the four Markov chain models eventually converge to a more accurate prediction. Finally, in Mahama, it is not immediately clear which models yield more accurate results, although the Graph-Adjusted and Time-Adjusted Models more consistently produce smaller ARDs.

This camp specific ARD is the main tool used in [8] for quantifying the success of their model of the Burundi crisis and thus 7 allows for a direct comparison with the results of [8]. Moreover, an informal inspection of the ensemble of camp ARD graphs indicates similar conclusions as found in the main text, specifically that the Markov chain models provide a modest improvement in the fit to the data compared to Flee and that the Markov Chain: Camp Adjusted model seems the best fit out of all of the models studied. In the main text these conclusions were quantitively supported by an analysis of the total ARD and Averaged ARGs, see Figures 3 & 6 and Tables 1 & 2.