Commuter Count:
Inferring Travel Patterns from Location Data

Nathan Musoke1,2, Emily Kendall1, Mateja Gosenca1,3,
Lillian Guo1, Lerh Feng Low1, Angela Xue1,
and Richard Easther1

1Department of Physics, University of Auckland, New Zealand

2Department of Physics and Astronomy, University of New Hampshire, USA

3Faculty of Physics, University of Vienna, Boltzmanngasse 5, 1090 Vienna,
Austria

For correspondence<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
##  1 Introduction

The movement of people between geographical regions and their personal
interactions are key determinants of the spread of pathogens such as Covid-19.
While interpersonal connections occur on scales ranging from individual
households to international travel, interactions between people in the course
of their daily routines provide a key “meso” layer in any detailed analysis of
pathogenic transmission.

The accumulation and analysis of data on the daily activities of individuals
has privacy implications and commercial sensitivities, creating (entirely
legitimate) barriers to its use by modelers. However, while it is unlikely
that detailed trajectories of individuals through the course of a day will be
shared outside of tightly controlled environments, aggregated spatio-temporal
data can often be made available.

In this Working Paper we analyse strategies for using aggregated spatio-
temporal population data acquired from telecommunications networks to infer
travel and movement patterns within regions. Specifically, we focus on hour-
by-hour cellphone counts for the SA-2 geographical regions covering the whole
of New Zealand [12] and base our work on algorithms described by Akagi et al.
[1]. This Working Paper describes the implementation of these algorithms,
their ability to yield inferences based on this data to build a model of
travel patterns during the day, and lays out opportunities for future
development.

Our testing data set consists of cellphone counts during January and February
2019 and 2020, where counts are given for individual New Zealand SA-2
geographical regions on an hourly basis. For reference, there are 2253 SA-2
regions in New Zealand. The regions vary in area, such that their typical
population is on the order of a few thousand people. The Greater Auckland
region contains approximately 600 SA-2s whereas in remote parts of the South
Island a similar geographical area might correspond to just handful of
SA-2s.111There are several exceptions, including offshore islands, which are
very thinly populated. This approach also implicitly assumes that cellphone
counts are a good proxy for the locations of the majority of the population.

We focus on the two algorithms, developed by Akagi and colleagues, referred to
as the ‘exact’ and ‘approximate’ methods. These algorithms use hour-by-hour
population counts to estimate bidirectional flows between pairs of
geographical regions. Long-distance travel over short time periods is
penalised by a simple function of the physical separation between regions.
Furthermore, a strict upper bound can be applied to the distance plausibly
travelled between successive time steps, so that not all possible region
pairings are viable. The algorithms adapt naturally to “real” geographies with
complex shapes (rather than a grid-based geometry) and data in which the total
number of individuals is not constant, due to phones being turned on or off or
moving in and out of coverage areas.

The motivation for this work was to facilitate analyses of the spread of
Covid-19 in New Zealand. However, the treatment here can be applied to any
number of tasks requiring models of population flows. This work investigates
these algorithms and extends our understanding of their properties and
limitations.

Having implemented both the exact and approximate algorithms, we test the
consistency of their outputs and find limitations and sensitivities to input
parameters which are common to both algorithms. We also identify and address
issues that arise when the number of people leaving a given region is roughly
similar to the number of destinations available to them, so that the expected
number of people moving between many pairs is less than one. In addition we
have developed a simulator which generates synthetic data that allows the
algorithms to be explored without access to cellphone counts and the
underlying individual trajectories, facilitating additional verification
strategies.

Our implementation of the exact algorithm is computationally efficient; far
more so than originally reported. In particular, we can “solve” for the
Greater Auckland and Waikato regions (encompassing the three Auckland area
District Health catchments) in tens of seconds of walltime on a laptop.

This Working Paper is structured as follows. Section 2 provides a quick
quantitative survey of the cellphone count data utilised in the inference
process. Section 3 discusses the construction of a likelihood function
characterising the probabilities of transitions between regions and Section 4
summarises the exact Akagi algorithm to maximise this likelihood, describes
its implementation in a computationally efficient Python code222Available at
https://github.com/auckland-cosmo/FlowStock, and identifies issues with its
application to this problem. Section 5 introduces a simple data simulator used
to validate the code, while Section 6 looks at the sensitivity of the results
to key parameters in the model. Section 7 describes the alternative
approximate algorithm proposed by Akagi et al. and contrasts its output to the
exact algorithm. In Section 8 we sketch an approach to compare our results to
external sources of commuter information. Finally, Section 9 provides a brief
summary of our experiences and identifies areas for further exploration. We
discuss possible improvements and extensions to the algorithms, and highlight
issues with these algorithms that might limit their utility.

##  2 Cellphone Count Data

Figure 1: Representative counts throughout a month. There are clearly large
daily commutes in and out of the central-city SA-2 region Auckland-University,
anti-correlated with flows to the residential area Balmoral. There is a
discernible difference between workday and weekend counts. Inspecting data
from Puketona-Waitangi, containing the Waitangi Treaty grounds, one can see a
significant increase in the lead up to February 6th.  Figure 2: Hourly
differences in the count in the Auckland-University area during the 26th to
28th of February 2020. One can see a sharp morning rush hour and less
pronounced evening rush hour. There is an anomaly at midnight on the 28th.
Such features are common at midnight and appear to be artifacts associated
with the capture and processing of the data by the telecommunications
providers.

Our analysis uses aggregated cellphone count data gathered from New Zealand
telecommunications companies. In particular, this data gives hourly counts of
the number of cellphones active in each SA-2 region. Note that the term
‘active’ applies to all cell phones which are powered on and detectable by the
cell network; a cell phone does not need to be in use to be considered active.
Within this data, it is possible to clearly discern patterns in population
flow, for example during weekends, holidays, or large gatherings. Figure 1
provides some representative examples.333February 6th is a public holiday in
New Zealand, during which there is often a large gathering at Puketona-
Waitangi to commemorate the signing of the Treaty of Waitangi.

It should be noted that each cell phone is counted in only one SA-2 region per
hour. This is reflected by the conservation of the total cell phone count over
time. Indeed, while a cell phone may be in range of multiple cell towers at
any given moment, it will only use a single tower for communication at any one
time, as determined by the relative signal strength. Hence, when the
instantaneous count is performed on an hourly basis, each cell phone is
associated with only one tower/SA-2 region. As the hourly data represents
instantaneous counts, movements between SA-2 regions/cell towers on timescales
smaller than one hour are not captured.

While most of the adult population is likely to carry a cellphone with them on
a day-to-day basis, there is no guarantee that cellphone counts map uniquely
to individuals or that the mapping is unbiased. Indeed, we expect that certain
demographics — e.g. the very young or very old, and the economically deprived
— may be missed in this data. Furthermore, populations with 0 or multiple
phones will be heavily biased in age, social class, and other areas that are
correlated with infection risk factors. Unfortunately, the currently available
data on cell phone access is not sufficiently detailed to incorporate into our
modelling at this time. While some relevant 2018 census data is available [8],
it only provides information on access to telecommunication systems at the
household, rather than individual, level. Furthermore, the census data
includes no information for 7.7% of households. While a detailed study of cell
phone ownership is outside of the scope of this work, it is expected that data
from future national surveys may improve our ability to correlate the
movements of cell phones with the movements of individual persons.

Finally, we also note that the data exhibits frequent discontinuities in
regional counts at midnight, when cell tower data is apparently “rebaselined”,
as shown in Figure 2. However, since our focus will be on movements during the
working day this is of limited relevance to our analysis.

##  3 Log-Likelihood and Likelihood Gradient

$V$ |  set of regions
---|---
$n$ |  number of regions, $|V|$
$T$ |  number of snapshots
$\mathbf{N}$ |  $n\times T$ matrix of counts in regions at each snapshot
$\mathbf{d}$ |  $n\times n$ matrix of distances $d_{ij}$ from region $i$ to $j$
$K$ |  distance cutoff
$V$ |  set of regions
$\Gamma_{i}$ |  set of neighbours of region $i$; $\\{j\in V|d_{ij}\leq K\\}$
${M}_{tij}$ |  the number of people who move from $i$ to $j$ between $t$ and $t+1$; $\mathbf{M}$ is a $(T-1)\times n\times n$ array
$\pi_{i}$ |  departure probability of region $i$
$s_{i}$ |  gathering scores for region $i$
$\beta$ |  scalar distance weighting
$\theta_{ij}$ |  probability for a person to move from region $i$ to $j$ between snapshots
$\mathcal{C}(\mathbf{M};\mathbf{N})$ |  cost function to enforce number conservation
$\lambda$ |  weighting of cost function
$\mathcal{L}(\mathbf{M},\bm{\pi},\mathbf{s},\beta;\mathbf{N},\lambda,\mathbf{d},K)$ |  likelihood of $\mathbf{M}$, $\bm{\pi}$, $\mathbf{s}$, and $\beta$ given data $\mathbf{N}$ and assumptions $\lambda$, $\mathbf{d}$, $K$
$\epsilon$ |  convergence threshold for iterative optimisation
Table 1: Symbols used in the text. Bold symbols are non-scalar quantities.

Following Akagi et al., we introduce a probability of transition between
different regions,

$P(\mathbf{M}|\mathbf{N},\bm{\theta})=\sum_{t=0}^{T-2}\sum_{i\in
V}\left(\frac{N_{ti}!}{\prod_{j\in\Gamma_{i}}{M}_{tij}!}\prod_{j\in\Gamma_{i}}\theta_{ij}^{{M}_{tij}}\right).$
(1)

Here $N_{ti}$ denotes the observed number of people in region $i$ at step $t$
(the algorithms can consider multiple time slices), which is provided as input
data. The number of transitions from region $i$ to $j$ at step $t$ is
represented by $M_{tij}$. The $M_{tij}$ are the quantities we seek to
estimate444Note that $T$ represents the total number of time slices, such that
there are $T-1$ time steps between the slices, labelled from $t=0$ to
$t=T-2$.. For each starting region $i$, the set of possible destination
regions is denoted by

$\Gamma_{i}=\\{j\in V|d_{ij}\leq K\\}$ (2)

where $d_{ij}$ is the distance from region $i$ to region $j$; $K$ is a cutoff
distance beyond which people are assumed not to travel in a single time
step.555We assume that the distance metric $d$ corresponds to the centroid-to-
centroid distance between SA-2 regions. Centroid coordinates are available at
https://datafinder.stats.govt.nz/layer/93620-statistical-area-2-2018-centroid-
true/. The probability of a person in region $i$ at time $t$ moving to region
$j$ at time $t+1$ is then $\theta_{ij}$. In general, this probability will be
dependent on the time of day. For example, commuter traffic into and out of
central business districts tends to reverse from morning to evening. It is
therefore important that the estimation algorithm be applied across time
periods in which transition probabilities may be assumed to be roughly
constant.

The algorithm requires an assumption for the transition probability, which is
taken to be

$\theta_{ij}=\begin{dcases}1-\pi_{i}&\text{\ if }i=j\\\
\pi_{i}\dfrac{s_{j}\exp(-\beta
d_{ij})}{\sum_{k\in\Gamma_{i}\setminus\\{i\\}}s_{k}\exp(-\beta
d_{ik})}&\text{\ if }i\neq j\end{dcases}\,,$ (3)

where the $\pi_{i}$ are components of the vector $\bm{\pi}$ of length $n$
which describes the probability of a person leaving their current region.
Their possible destinations are weighted by another $n$-vector, $\mathbf{s}$,
where $s_{j}$ describes the tendency for people to gather in region $j$. For
example, regions within the central business district would be expected to
have a strong tendency to attract commuters during the morning rush hour.
Following Akagi et al., we include an exponential penalty on long-distance
travel, but note that other forms of penalty are possible.666As we see below,
$\beta$ is one of the parameters we adjust to optimise the fit. In some cases
the optimal value of $\beta$ was negative, but often for unrealistically small
regions — and there are also more possible pairings at greater distances. We
experimented with the choice $e^{-\beta_{1}d_{ij}+\beta_{2}d_{ij}^{2}}$ but
did not pursue it in detail. Finally, note that $\mathbf{s}$ has an arbitrary
overall normalisation.

Akagi et al. obtain a log-likelihood function from Equation 1,

${\cal{L}}^{\prime}=\mathcal{L}_{0}+\mathcal{L}_{1}+\mathcal{L}_{2}-\frac{\lambda}{2}{\cal{C}({\mathbf{M},\mathbf{N}})},$
(4)

where the individual components are given by:

$\displaystyle\mathcal{L}_{0}=\sum_{t=0}^{T-2}\sum_{i}\log(1-\pi_{i}){M}_{tii},$
(5)
$\displaystyle\mathcal{L}_{1}=\sum_{t}\sum_{i}\sum_{j\in\Gamma_{i}\backslash\\{i\\}}\left(\log(\pi_{i})+\log(s_{j})-\beta
d_{ij}-\log\sum_{k\in\Gamma_{i}\backslash\\{i\\}}s_{k}e^{-\beta
d_{ik}}\right){M}_{tij},$ (6)
$\displaystyle\mathcal{L}_{2}=\sum_{t}\sum_{i}\sum_{j\in\Gamma_{i}}(1-\log{M}_{tij}){M}_{tij},$
(7)
$\displaystyle{\cal{C}({\mathbf{M},\mathbf{N}})}=\sum_{t=0}^{T-2}\sum_{i}{\left(N_{ti}-\sum_{j}{M}_{tij}\right)}^{2}+\sum_{t=0}^{T-2}\sum_{i}{\left(N_{t+1,i}-\sum_{j}{M}_{tji}\right)}^{2}.$
(8)

Stirling’s approximation for factorials is used in the course of this
derivation; we will revisit this choice in Section 6.1. The diagonal component
of the $t$-th transition matrix $M_{tii}$ corresponds to the population that
does not leave block $i$ at step $t$. The cost function
$\mathcal{C}(\mathbf{M},\mathbf{N})$ is a soft enforcement of number
conservation and this is the only place where the overall size of the
population enters the likelihood, rather than dimensionless transition rates.
The strength of the cost function is controlled by the parameter $\lambda$.

We estimate flows by maximizing $\mathcal{L}$ with respect to the $n^{2}$
components of $\mathbf{M}$ (per time step), the $n$ components of $\bm{\pi}$
and $\mathbf{s}$, and the scalar $\beta$. The distance cutoff can fix some
components of $\mathbf{M}$ to zero but this will not always result in a
meaningful simplification of the optimisation problem. For instance, the
Auckland region directly couples in excess of 500 SA-2 blocks, creating
approximately 250,000 pairs, each of which has a corresponding $M_{tij}$.
Consequently, the application of this algorithm to a realistic problem
involves estimating values for $10^{5}$ to $10^{6}$ variables.

We perform the optimisation with the SciPy [2, 14] implementation of the
L-BFGS-B algorithm. By default, derivatives of the target function are
evaluated via differencing, requiring multiple evaluations of the likelihood.
Since the complexity of the likelihood and the number of free parameters both
grow with the number of possible pairs the optimisation quickly becomes
numerically challenging. However, we can greatly improve performance by
supplying analytic derivatives as the ${M}_{tij}$ do not appear in complicated
combinations within the likelihood.

After some calculation, we find that the derivatives of the terms in Equation
4 are

$\displaystyle\frac{\partial\mathcal{L}_{0}}{\partial{M}_{tij}}=\begin{cases}\log(1-\pi_{i})&i=j\\\
0&i\neq j\end{cases}$ (9)
$\displaystyle\frac{\partial\mathcal{L}_{1}}{\partial{M}_{tij}}=\begin{cases}0&i=j\\\
\log(\pi_{i})+\log(s_{j})-\beta d_{ij}-\log\sum_{k\in\Gamma_{i}\backslash
i}s_{k}e^{-\beta d_{ik}}&j\in\Gamma_{i}\setminus\\{i\\}\\\
0&j\notin\Gamma_{i}\end{cases}$ (10)
$\displaystyle\frac{\partial\mathcal{L}_{2}}{\partial{M}_{tij}}=\begin{cases}-\log{M}_{tij}&j\in\Gamma_{i}\\\
0&j\notin\Gamma_{i}\end{cases}$ (11)
$\displaystyle\frac{\partial\mathcal{C}}{\partial{M}_{tij}}=-2\left(N_{ti}-\sum_{l}M_{til}\right)-2\left(N_{t+1,j}-\sum_{l}M_{tlj}\right)$
(12)

While computing the likelihood requires summing $\mathcal{O}(n^{2})$ terms,
the derivative of the cost function requires summing $\mathcal{O}(n)$ terms,
and each of other terms in the derivative is a single term from the sums in
$\mathcal{L}$. Consequently, evaluating the derivative of $\mathcal{L}$ with
respect to each of the $n^{2}$ components of $\mathbf{M}$ will involve
$\mathcal{O}(n\times n^{2})$ operations whereas approximating them numerically
would involve $n^{2}$ evaluations of the likelihood, for a total cost of
$\mathcal{O}(n^{2}\times n^{2})$.

We may further improve computational efficiency in evaluating both
$\mathcal{L}$ and $\partial\mathcal{L}/\partial{M}_{tij}$: when optimising
with respect to $\mathbf{M}$, the $\log$ in Equation 5 and bracketed term in
Equation 6 do not change, and can be precomputed, offering a significant
improvement in efficiency over millions of evaluations.

##  4 “Exact” maximisation algorithm

The “exact” maximisation algorithm described by Akagi et al. requires an
iterative maximisation, looping over three separate maximisations until the
relative difference in $\mathcal{L}$ changes by less than $\epsilon$, an
adjustable parameter. We have implemented the exact algorithm as follows:

  1. 1.

Initialise $\mathbf{M}$, $\bm{\pi}$, $\mathbf{s}$, $\beta$.

This step is unspecified in [1], but the way it is done has a significant
impact on the results of the algorithm. We discuss this further in Section
6.2.

  2. 2.

Loop over steps below until the relative difference in $\mathcal{L}$ changes
by less than $\epsilon$.

    1. (a)

Maximise $\mathcal{L}$ with respect to $\mathbf{M}$ while keeping $\bm{\pi}$,
$\mathbf{s}$ and $\beta$ constant.

    2. (b)

Maximise $\mathcal{L}$ with respect to $\bm{\pi}$ while keeping $\mathbf{M}$,
$\mathbf{s}$, $\beta$ constant, via the exact expression from Akagi et al.

$\pi_{i}=\frac{\sum_{t}\sum_{j\in\Gamma_{i}\setminus\\{i\\}}{M}_{tij}}{\sum_{t}\sum_{j\in\Gamma_{i}}{M}_{tij}}\,.$
(13)

    3. (c)

Iteratively optimise $\mathcal{L}$ with respect to $\mathbf{s}$ and $\beta$.

In contrast with Akagai et al., who use the Minorisation-Maximisation
algorithm for this step, we optimise the $\mathbf{s}$ and $\beta$ dependent
part of $\mathcal{L}$ directly. Our experience was that the former approach
can become “stuck” during an evaluation.

The only part of $\mathcal{L}$ that depends on $\mathbf{s}$ and $\beta$ is
$\mathcal{L}_{1}$ and it can be rearranged into a target function $f$ defined
by

$\displaystyle f=\sum_{i\in
V}\left(A_{i}\log(s_{i})-B_{i}\log\left(\sum_{k\in\Gamma_{i}\setminus\\{i\\}}s_{k}\exp(-\beta
d_{ik})\right)\right)-\beta D$ (14) $\displaystyle
A_{i}=\sum_{t}\sum_{j\in\Gamma_{i}\setminus\\{i\\}}{M}_{tji}$ (15)
$\displaystyle B_{i}=\sum_{t}\sum_{j\in\Gamma_{i}\setminus\\{i\\}}{M}_{tij}$
(16) $\displaystyle D=\sum_{t}\sum_{i\in
V}\sum_{j\in\Gamma_{i}\setminus\\{i\\}}d_{ij}{M}_{tij}\,.$ (17)

The derivation of $A_{i}$ requires reordering the sum containing $\mathbf{s}$.
This resummation obscures the scale-independence of $\mathbf{s}$ seen in
Equations 1 and 4, and is only valid when the matrix $\mathbf{d}$ of distances
$d_{ij}$ is symmetric.777The $d_{ij}$ is effectively a cost function for
travel between Block-$i$ and Block-$j$. We have assumed this is symmetrical
and $d_{ij}=d_{ji}$ but in principle this could be (for example) time-
dependent and incorporate congestion related delays. We do not consider these
possible asymmetries in $d$. We proceed as follows:

      1. i.

Optimise $f$ with respect to $\mathbf{s}$. There is a closed form for this:

$s_{i}=\frac{A_{i}}{\sum_{k}C_{k}\exp(-\beta d_{k_{i}})}$ (18)

      2. ii.

Normalise $\mathbf{s}$ with

$\mathbf{s}\mapsto\frac{\mathbf{s}}{\max(\mathbf{s})}\,.$ (19)

This is done to avoid numerical problems where $|\mathbf{s}|\to 0$ otherwise.

      3. iii.

Maximise $f$ with respect to $\beta$. This maximisation is done with the
bounded Brent algorithm.

We found that this optimisation of $\mathbf{s}$ and $\beta$ would occasionally
enter a closed loop. When this happens, we terminate the optimisation of
$\mathbf{s}$ and $\beta$ and return to the optimisation of $\mathbf{M}$ and
$\bm{\pi}$ before trying again.

We note the similarity of the procedure described here to the well-known
Expectation Maximisation (EM) algorithm [3]. The EM algorithm is a method for
performing maximum likelihood estimation in the presence of latent variables
and has broad applicability to diverse fields from computational biology to
machine learning [5, 4]. The EM algorithm works by iteratively improving
parameter value estimates through alternating between an expectation step and
a maximisation step. In the expectation step, a function for the expectation
value of the log-likelihood is computed using the existing estimations for the
parameter values. In the subsequent maximisation step, new parameter values
are computed which maximise the value of the previously determined function.
This process is then iterated until the desired convergence criteria are met.
An adaptation of the EM algorithm which closely resembles our approach is
known as Expectation Conditional Maximisation (ECM) [7]. In this case, each
maximisation step is subdivided into a series of conditional maximization
steps in which maximisation with respect to each parameter is undertaken
individually, while all other parameters remain fixed. A detailed comparison
between the efficacy of the algorithm implemented in this work and other
variants of EM/ECM is out of scope here, but warrants further investigation
going forward.

The majority of the implementation of this “exact” maximisation algorithm is
in idiomatic Python and Numpy [13]. Some calculations make use of Numba [6],
but ultimately this was not a major performance gain over vanilla Numpy.
Including the analytic functions in Equations 9, 10, 11 and 12 for the
derivative of the likelihood improved performance by multiple orders of
magnitude.

Figure 3 shows wall clock computational time for a representative test problem
based on synthetic data. While a plateau is observed in Figure 3 when the
number of regions exceeds $\sim 120$, we would of course expect further
increases in run time for much larger data sets. The details of this will
depend upon multiple factors such as NumPy’s usage of multiple CPUs for
certain operations, and the available memory. In practice, estimating
movements within the approximately $800$ SA2 blocks in the Auckland and
Waikato regions took $\sim 30$ seconds on a laptop; this is consistent with
synthetic data. Consequently, the numerical performance of our implementation
of the exact algorithm presents no significant challenges for any currently
envisaged applications, and appears to improve significantly on that reported
by Akagi et al., presumably thanks to the use of the analytic derivatives.

Figure 3: Plot of the run time against number of regions. The shaded bands
represent the standard deviation across 18 runs with 3 random seeds and 2
noise amplitudes for the synthetic data, and 3 choices of $\lambda$ in the
solver. Clearly, demanding higher precision increases the run time but it
remains manageable even as the number of regions grows beyond 100. All
simulations were run on a consumer grade computer.

##  5 Synthetic data

We do not have the true values of $\mathbf{M}$ for the cellphone data, whereas
Akagi et al. had access to trajectory information. However, we can test
against synthetic data that strictly matches the assumptions laid out in
Section 3.

We specify the number of regions we want to simulate, the distances $d_{ij}$
between them, cutoff distance $K$ and distance weighting $\beta$. Then we
stipulate vectors of gathering scores $\mathbf{s}$ and departure probabilities
$\bm{\pi}$ corresponding to each region. From this, the simulator calculates
the set of possible destinations $\Gamma_{i}$ of each region and probabilities
$\theta_{ij}$ for moving to each from Equations 2 and 3.

We specify an initial distribution of people in each region as a vector
$N_{0}$ and number of time steps to take. Then for each time step $t$ and
region $i$, the simulator loops over each of the people currently in region
$i$ and assigns them to be in region $j$ at time $t+1$ with probability
$\theta_{ij}$. This defines a “true” ${M}_{tij}$ Optionally, the simulator
also randomly adds or removes people before calculating where they go. This
allows us to test against scenarios that do not conform exactly to the
assumptions in Section 3.

Figure 4: From left to right: true $s$, true $\pi$, initial counts, and final
counts for simulated synthetic data. This data has 9 regions, each with an
initial count of $N_{0}=1,000,000$. People leave each region with probability
$\pi$ and each region has a gathering score $s$. One can see that the region
with higher $\pi$ has a net loss in $N_{1}$. The regions with larger $s$ have
correspondingly larger net gains.

Figure 4 shows a deliberately simple setup: an initially uniform population
distributed among 9 regions arranged in a regular $3\times 3$ grid. The
distance between the centres of the outermost cells on each side of the grid
is therefore 2 units.

We set the distance threshold $K=2$ and $\beta=1$. Because the grid (expressed
in terms of the centroids) has a width of 2, only corner to corner travel is
disallowed with $d=2\sqrt{2}>K$. The departure probabilities $\bm{\pi}$ are
sampled uniformly from $0.01$–$0.02$, other than the central region which has
a probability of $0.1$. This higher departure probability is evident in the
populations after one time step; the central region has a larger net loss of
people than the outer regions. The gathering scores $\mathbf{s}$ are set to 1
other than 4 regions shown in the left panel of Figure 4. The gathering scores
have the expected effect; regions with larger $s$ have correspondingly larger
net gains.

Figure 5: Scatter plots of the true and estimated values of $s_{i}$, $\pi_{i}$
and ${M}_{tij}$. The accuracy of the $s$ and $\pi$ estimates are relatively
good, but there are a number of the transition matrix elements $M$ that are
severely underestimated. The large elements of ${M}_{tij}\sim 10^{6}$ are in
the diagonals; these are people who did not move. Figures such as this are
presented as 2-D histograms in the following sections, where there are too
many points for a sensible scatter plot.

Figure 5 shows the results of applying our implementation of the exact
algorithm to the data described above. It is able to get good estimates of the
gathering scores and departure probabilities. However, the estimates of
${M}_{tij}$ and $\beta$ are poor. There are a number of transitions that are
severely underestimated. In addition, $\beta$ is estimated at 0.08 rather than
the true value of 1.0. Poor estimates of $\beta$ are a recurring theme in the
following sections.

##  6 Implementation and Validation

We now consider the stability of the results against changes in free
parameters in the algorithm, and specific issues that arise when we apply the
Akagi et al. algorithms to the present use-case.

###  6.1 Scaling

As mentioned in Section 3, Equation 4 assumes that Stirling’s approximation
$\log{n!}\approx n\log n-n$ applies to the elements of the transition matrix,
or ${M}_{tij}\gtrsim 1$. However, this assumption is violated by fairly
typical real-world data. SA2 regions generally contain $\mathcal{O}(1,000)$
people and if $\mathcal{O}(100)$ people enter or leave in an hour with
$\mathcal{O}(100)$ allowed destinations and origins, some transitions will
necessarily involve “fractional” numbers of people. These fractional numbers
of people should be interpreted as probabilities.

Figure 6: Histograms comparing computed $M_{tij}$ for different population
scalings. The data is for transitions between the 798 unique SA-2s in the
regional councils of Auckland Region and Waikato. The plots show the counts of
$M_{tij}$ for pairs of scalings; with perfect agreement all elements would lie
on the diagonal, up to the scatter arising from the large number of near-
degenerate solutions to the maximisation problem. The top panel compares the
raw counts (a scaling of 1) with a scaling of 1000 (y-axis). The bottom panel
compares a scaling of 1000 (x-axis) and 10,000 (y-axis).

We have found that the inapplicability of Stirling’s approximation can be
ameliorated by scaling up the number of people to the point where all allowed
transitions have $M_{tij}\gtrsim 1$. The cost function, Equation 8 is
quadratic in this scale factor, so one must simultaneously rescale $\lambda$
to compensate. For sufficiently large scalings the results become scaling
independent. We checked that this is true by comparing the results for the
SA2s contained in the combined Auckland Region and Waikato on February 18th
between 7am and 9am. There are 798 regions in total, but the small number with
less than 100 people are dropped from the analysis of scaling by 1000 and
10,000, as shown in Figure 6. This strategy is not necessarily perfect — the
${M}_{tij}\log({M}_{tij})$ term in $\mathcal{L}_{2}$ is non-linear in
${M}_{tij}$ and requires more detailed analysis — but will be more robust than
using unscaled populations. All other results shown in the following sections
use a scaling large enough to ensure that the computed ${M}_{tij}\gtrsim 1$
and these are then rescaled back to their original values.

###  6.2 Repeatability and Initial conditions

The solver needs initial conditions that do not represent pathological
states.888Note that ‘initial conditions’ does not refer to values at $t=0$ but
to the initial guesses for $\beta$, $\mathbf{s}$, $\bm{\pi}$, and the entire
$\mathbf{M}$ matrix, prior to optimisation. As an initial guess we make the
following “static” choice

$\displaystyle\pi_{i}=0.02$ (20) $\displaystyle s_{i}=0.02$ (21)
$\displaystyle\beta=\frac{50}{\max(d)}$ (22) $\displaystyle
M_{tij}=\begin{cases}N_{ti},&\text{for $i=j$ }\\\ 0,&\text{for $i\neq
j$}\end{cases}\,,$ (23)

where the last line implies that no-one moves between blocks. This is not
self-consistent, since if the off-diagonal $M_{tij}=0$, the $\pi_{i}$ and
$s_{i}$ should also vanish, but these initial conditions can yield plausible
results. Varying the starting position causes the algorithm to converge on
different maxima.

We first test the sensitivity to initial conditions by adding a random scatter
to the initial guess:

$\displaystyle M_{tij}=\begin{cases}N_{ti}+\delta_{tii},&\text{for $i=j$ }\\\
\delta_{tij},&\text{for $i\in\Gamma_{i}\setminus\\{i\\}$ }\\\ 0,&\text{for
$i\notin\Gamma_{i}$}\end{cases}$ (24)

where $\delta_{tij}$ is sampled uniformly from the range $[0,N_{ti})$. In
Figure 7 we show that this scatter in the initial conditions does not have a
drastic impact on the output by analysing data from SA2s contained in the
combined Auckland Region and Waikato regions, on February 18th at 7am 8am, and
9am. We quantify this sensitivity by computing the mean and standard deviation
for the values of each of the ${M}_{tij}$ from a sequence of 20 runs with
different random perturbations to the original initial condition Equation 23.
We find that the ratio of the standard deviation to the mean is small for the
vast majority of ${M}_{tij}$ for the cases we consider.

Figure 7: Top: Histogram of the normalised standard deviation
$\mathrm{std}({M}_{tij})/\bar{M}_{tij}$ of $M$ values. The mean is with
respect to 20 runs, each with the initial conditions of Equation 24 and
different seeds for the random jitter; in most cases
$\mathrm{std}(M_{tij})/\bar{M}_{tij}$ is significantly less than unity.
Bottom: Two-dimensional histogram of the same data. When the standard
deviation is below the red line it is less than the mean value; $M_{tij}$
above this line have large scatter between runs. Note the logarithmic colour
scale on this plot; the most common points and the largest $M$ values are
below the line. Data is for SA2s contained in the combined Auckland Region and
Waikato, on February 18th at 7am 8am, and 9am.

We also consider a “moving” initial guess,

$M_{tij}=\begin{cases}N_{ti},&\text{for $i=j$ }\\\
\frac{\left|N_{ti}-N_{t+1,i}\right|}{|\Gamma_{i}\setminus\\{i\\}|},&\text{for
$j\in\Gamma_{i}\setminus\\{i\\}$}\\\ 0,&\text{for
$i\notin\Gamma_{i}$}\end{cases}\,.$ (25)

This encodes an expectation that most people stay where they are and that the
number of people moving out of a region is on the order of the change in its
population (regardless of whether that change is a net inflow or outflow).

In Figure 8 we compare the two initial conditions choice described above. We
use data from the 63 most populated areas in Southland Region, on 11 February
2020 at 6am, 7am, 8am, 9am and 10am. There is a clear discrepancy when
$\epsilon=10^{-2}$ but moving to a more stringent $\epsilon=10^{-4}$
eliminates much of this bias.

Figure 8: Histograms of inferred $\mathbf{M}$ values with different initial
conditions. The top panel has $\epsilon=10^{-2}$ and the bottom has
$\epsilon=10^{-4}$. Static initial conditions ($x$-axis) start with diagonal
transition matrices; i.e. no-one moves, as in Equation 24; the $y$-axis have
initial conditions for which many people move, as in Equation 25. With a loose
convergence parameter the final result reflects the initial choice; setting
$\epsilon=10^{-4}$ eliminates most sensitivity to initial conditions. Data is
from the 63 most populated areas in Southland Region, on 11 February 2020 at
6am, 7am, 8am, 9am and 10am.

###  6.3 Sensitivity to $\epsilon$ and $\lambda$

The values $\epsilon$ and $\lambda$ have an impact on the output of the
algorithm. We used the normalised absolute error (NAE) to quantify the error
in these fits, which is defined to be

$\frac{\sum\limits_{t,i,j}\left|M^{*}_{tij}-M_{tij}\right|}{\sum\limits_{t,i,j}M^{*}_{tij}},$
(26)

where $M^{*}_{tij}$ represents the ‘true’ $M$ values from the simulated data.
We note that the NAE may be misleading in cases where there are a small number
of regions for which the $M_{tij}$ are highly inaccurate if these regions have
relatively large populations.

We examined the impact of $\epsilon$ and $\lambda$ by running the exact
estimator on simulated data, as in Section 5. We assume $15^{2}$ cells
distributed on a regular $2\times 2$ grid and a distance cutoff $K=1.5$, so
that each region has $100$ to $225$ possible destinations. The initial number
of people, gathering scores and leaving probabilities in cell $i$ were set to

$\displaystyle
N_{0,i}=\nu\exp\left(-{(\sqrt{x_{i}^{2}+y_{i}^{2}}-r_{0})}^{2}\right)$ (27)
$\displaystyle s_{i}=\exp(-4(x_{i}^{2}+y_{i}^{2}))$ (28)
$\displaystyle\pi_{i}=\frac{1}{10}\frac{N_{0,i}}{\max_{j}(N_{0j})}$ (29)
$\displaystyle\beta=1$ (30)

where $r_{0}=0.8$. The gathering score is high at the center, and the
departure probability is proportional to the initial number of people in a
cell. There is a higher density of people in a ring of radius $r_{0}$ around
the center. This is intended to be roughly analogous to people migrating from
the outskirts of a city into the center. We allowed a 10% error in the number
of people at each location during each time step.

The results are shown in Figure 9. The absolute variance in the NAE as a
function of both $\epsilon$ and $\lambda$ is not large. Counterintuitively, we
found that smaller values of $\epsilon$ do not necessarily give more accurate
results by this measure, but the differences are not significant. There is no
obvious choice of $\lambda$; large values heavily penalise solutions where
summing the people going in and out of regions does not match the known data
$N$. There are also numerical issues introduced by large $\lambda$; these seem
much like the issues introduced with very small $\epsilon$. Small values of
$\lambda$ allow proposed solutions to have large violations of number
conservation. Given that the real-world data is known to have imperfect number
conservation, some deviation should be allowed and a middle ground should be
found. The bottom panel in Figure 9 confirms this intuition.

Figure 9: The value of the NAE as compared to simulated data estimator
parameters $\epsilon$ (top) and $\lambda$ (bottom). The error bands come from
aggregating over 8 instances of $N$ and $\lambda$ (top) and $\epsilon$
(bottom). The blue solid line is the error when the simulated data conforms
exactly to the assumptions of the likelihood Section 3. The orange dashed line
assumes that there is an error of up to 10% at each step. Interestingly, the
estimator performs better on noisy data. Smaller values of $\epsilon$ do not
have a clear advantage when there is noise in the data. On the other hand,
$\lambda=10$ is better than 1 or 100.

##  7 Alternative Algorithm: Approximate Inference Method

Akagi et al. present an alternative method, which is billed as being less
computationally expensive. This concern is less pressing, given the speedup
applied to the exact algorithm. We have implemented a variation of this
algorithm in Python, with a few key differences. In particular, Akagi et al.
bin regions by their separations but we cannot adopt this approach, given the
irregular spacing of the SA-2 centroids.

###  7.1 Summary of alternative algorithm

We begin by defining the following parameters:

$X_{tij}\equiv M_{tji},\quad Y_{ti}\equiv\sum_{j\neq i}M_{tij},\quad
Z_{ti}\equiv M_{tii}.$ (31)

Using these parameters, $\bm{\pi}$ and $f(\mathbf{s},\beta)$ are given by:

$\pi_{i}=\frac{\sum\limits_{t}Y_{ti}}{\sum\limits_{t}(Y_{ti}+Z_{ti})},$ (32)

$f(\mathbf{s},\beta)=\sum_{t,i,j}(X_{tij}\log s_{i}-\beta
d_{ij}X_{tij})-\sum_{t,i}Y_{ti}\log\Big{(}\sum_{k\neq i}s_{k}\exp(-\beta
d_{ij})\Big{)},$ (33)

where $d_{ij}$ is the distance between centroids of SA-2 regions. We also
define parameters $\theta_{ij}$ and $\mu_{ij}$ as follows:

$\theta_{ij}=\begin{cases}1-\pi_{i},&(i=j)\\\
\pi_{i}\left(\frac{s_{j}\exp(-\beta d_{ij})}{\sum\limits_{k\neq
i}s_{k}\exp(-\beta d_{ij})}\right),&(i\neq j)\end{cases}$ (34)

$\mu_{ij}=\sum_{t}N_{tj}\theta_{ji}.$ (35)

Following Akagi et al., the approximate log likelihood is given by

$\displaystyle\mathcal{L}_{\text{approx}}=$
$\displaystyle\sum_{t,i,j}\big{(}X_{tij}\log(\mu_{ij})+X_{tij}-X_{tij}\log(X_{tij})\big{)}$
$\displaystyle+\sum_{t,i}\big{(}Y_{ti}\log(N_{ti}\pi_{i})+Y_{ti}-Y_{ti}\log(Y_{ti})\big{)}$
$\displaystyle+\sum_{t,i}\big{(}Z_{ti}\log(N_{ti}(1-\pi_{i}))+Z_{ti}-Z_{ti}\log(Z_{ti})\big{)},$
(36)

with the associated constraint function,

$C(X,Y,Z)=\sum_{t,i}\Big{(}|N_{ti}-(Y_{ti}+Z_{ti})|^{2}+|N_{t+1,i}-\sum_{j}X_{tij}|^{2}\Big{)}.$
(37)

We then have a log likelihood function for the final calculation of $M$ as
follows:

$\displaystyle\mathcal{L}_{\text{final}}=$
$\displaystyle\sum_{t,i}\log(1-\pi_{i})M_{tii}+\sum_{t,i,j}\big{(}M_{tij}-M_{tij}\log(M_{tij})\big{)}$
$\displaystyle+\sum_{t,i,j\neq i}\bigg{(}\log(\pi_{i})+\log(s_{j})-\beta
d_{ij}-\log\Big{(}\sum_{k\neq i}s_{k}\exp(-\beta d_{ik})\Big{)}\bigg{)},$ (38)

with associated constraint function:

$C(M)=\sum_{t,i}\Big{(}|N_{ti}-\sum_{j}M_{tij}|^{2}+|N_{t+1,i}-\sum_{j}M_{tji}|^{2}\Big{)}.$
(39)

The inference proceeds as follows:

  1. 1.

Initialise parameters $\mathbf{M}$, $X$, $Y$, $Z$, $\bm{\pi}$, $\mathbf{s}$,
and $\beta$,

  2. 2.

Maximise $\mathcal{L}_{\text{approx}}$ \- $\frac{\lambda}{2}C(X,Y,Z)$,

  3. 3.

Update $\bm{\pi}$,

  4. 4.

Update $\mathbf{s}$ and $\beta$ by maximising $f(\mathbf{s},\beta)$,

  5. 5.

Repeat 1 - 4 until specified convergence criterion is reached for the value of
the approximate log likelihood,

  6. 6.

Calculation of $\mathbf{M}$ through Maximising $\mathcal{L}_{\text{final}}$ \-
$\frac{\lambda}{2}C(\mathbf{M})$, using the final $\bm{\pi}$, $\mathbf{s}$,
and $\beta$ values calculated above.

When optimising $\mathcal{L}_{\text{approx}}$ and
$\mathcal{L}_{\text{final}}$, it is useful to define their analytic Jacobians,
as it is computationally expensive to compute approximate derivatives, along
with analytic Jacobians for the constraint. These are as follows:

$\displaystyle\frac{\partial\mathcal{L}_{\text{approx}}}{\partial X_{tij}}$
$\displaystyle=\log(\mu_{ij})-\log(X_{tij}),$ (40)
$\displaystyle\frac{\partial\mathcal{L}_{\text{approx}}}{\partial Y_{ti}}$
$\displaystyle=\log(N_{ti}\pi_{i})-\log(Y_{ti}),$ (41)
$\displaystyle\frac{\partial\mathcal{L}_{\text{approx}}}{\partial Z_{ti}}$
$\displaystyle=\log(N_{ti}(1-\pi_{i}))-\log(Z_{ti}),$ (42)

with constraint function derivatives:

$\displaystyle\frac{\partial\big{(}-\frac{\lambda}{2}C(X,Y,Z)\big{)}}{\partial
X_{tij}}$ $\displaystyle=\lambda\Big{(}N_{t+1,i}-\sum_{k}X_{tik}\Big{)},$ (43)
$\displaystyle\frac{\partial\big{(}-\frac{\lambda}{2}C(X,Y,Z)\big{)}}{\partial
Y_{ti}}$
$\displaystyle=\frac{\partial\big{(}-\frac{\lambda}{2}C(X,Y,Z)\big{)}}{\partial
Z_{ti}}=\lambda\Big{(}N_{ti}-(Y_{ti}+Z_{ti})\Big{)}.$ (44)

For the final log likelihood, we have:

$\displaystyle\frac{\partial\mathcal{L}_{\text{approx}}}{\partial M_{tii}}$
$\displaystyle=\log(1-\pi_{i})-\log(M_{tii}),$ (45)
$\displaystyle\frac{\partial\mathcal{L}_{\text{approx}}}{\partial M_{tij\neq
i}}$ $\displaystyle=\log(\pi_{i})+\log(s_{j})-\beta
d_{ij}-\log\Big{(}\sum_{k\neq i}s_{k}\exp(-\beta
d_{ik})\Big{)}-\log(M_{tij}),$ (46)

with constraint function derivatives:

$\frac{\partial\big{(}\frac{-\lambda}{2}C(M)\big{)}}{\partial
M_{tij}}=\lambda\Big{(}N_{ti}+N_{t+1,i}-\sum_{k}M_{tik}-\sum_{k}M_{tkj}\Big{)}.$
(47)

###  7.2 Performance of alternative algorithm

We implemented this algorithm in both Python 2 and Python 3, noting that the
former tends to outperform the latter, apparently due to bugs within the Numba
compiler in Python 3. Our implementation was first tested using synthetic
data. Using the NAE as a measure of the performance of the algorithm it was
found that for large data sets, it is beneficial to nest the main loop, as
described by steps 1 to 6 above, within an outer loop. This outer loop feeds
the calculated $\mathbf{M}$ values back as initial conditions in the
subsequent evaluation. For testing purposes, the outer loop is terminated
either when the NAE reaches a specified target value, or when successive loops
result in no further decrease in the NAE. This is only possible when the true
values are known, as in our test case. Hence, when applying this algorithm to
real-world data, one may choose to terminate the outer loop when the
successive change in $M_{tij}$ values reaches a certain threshold.

As an example, using simulated data with 225 regions over 3 time steps gave an
NAE of 0.046 after three iterations through the outer loop, compared to 0.154
with only one iteration. By comparison, the “exact” algorithm achieved an NAE
of 0.100, so that in this case the alternative algorithm appears to perform
better, though it does take more computation time. In this case we chose
$\lambda=10$, a convergence criterion of $0.001\%$ on the approximate log-
likelihood, and a tolerance $\texttt{ftol}=10^{-4}$ within
scipy.optimise.minimise for the $\mathbf{M}$ calculation.

We can also calculate the off-diagonal NAE, as the large values on the
diagonal can dominate the NAE, obscuring how well the algorithm is able to
identify regions with high gathering scores. In this case, the off-diagonal
NAE for the alternative algorithm was 0.279, compared with 0.558 for the
“exact” algorithm, again indicating more accurate reproduction of the input
data.999The synthetic data used in this test, along with the implementation
used to analyse it, are available in the code base within the folder ‘nae-
comparison’.

###  7.3 Discussion of alternative algorithm

Testing indicates that initialising the $\mathbf{M}$ arrays with the
corresponding $\mathbf{N}$ values on the diagonal and small random numbers on
the off-diagonal provides the best outcomes.101010The introduction of explicit
randomness in the $\mathbf{M}$ initialisation can make it difficult to compare
successive runs. To overcome this one may fix the seed of the random number
generator. The alternative inference algorithm runs through the entire
inference process multiple times, inputting the new $\mathbf{M}$ arrays as
initial conditions in each run. In some cases this leads to much improved
results, but can also result in an ‘overshoot’, whereby the off-diagonal
elements become too high.

The output is highly sensitive to the value of $\lambda$, which controls the
strength of the penalty terms. If $\lambda$ is too small, the algorithm tends
to overpopulate the off-diagonals. Conversely, if the value is too high, all
off-diagonal elements tend to zero. The optimal value of $\lambda$ varies on a
case-by case basis, making it difficult to guess a suitable value in advance.

In addition to the algorithm’s sensitivity to the $\mathbf{M}$ initialisation,
$\lambda$ value, and number of complete inference loops, one must also
consider the convergence criteria set for the approximate log-likelihood in
the inner loop, and tolerances set in the optimisation routines as well.
Tighter convergence constraints may increase computation time to an
unacceptable degree, or may preclude convergence entirely.

The original Akagi et al. treatment introduced this algorithm for its greater
efficiency, and it serves as a useful counterpoint to the “exact” version.
However, given its relative fragility with respect to control parameters and
the efficiency of our implementations it does not appear to offer a
significant, generic advantage.

##  8 Validation Test Case: Southland

We are unable to test the performance of the algorithms against real-world
trajectory information. However, one may gauge how well the algorithm captures
the essential features of daily travel by comparing its output to census data,
and we take a very cursory look at this problem here. We accessed the publicly
available self-declared commuting trends from the 2013 census using the Stats
NZ ‘Commuter View’ tool [9]. This tool presents the number of residents who
commute out of that region for work and the number that commute in. We can
then compare the trends in the census data to the sum of the off-diagonal
$\mathbf{M}$ matrix elements for outbound and inbound travel for each region
on a standard weekday, assuming most workers travel to work in a time period
from 6am to 10am.

For a simple test case, we have singled out the SA-2 regions which belong to
the wider Southland district, discarding regions with missing data at any of
the times of interest. Of the 65 SA-2 regions within the Southland district,
only 2 regions had incomplete data, both with populations of less than $10$.
This comparison is not exact. The subdivision of geographical regions within
the census data does not match the SA-2 regions used in the telecommunications
data so assumptions must be made when matching the data to our calculations.
Furthermore, this method does not capture the varying work hours of the
general populace, and is seven years out of date.

We ran the approximate inference algorithm for the telecommunications data
from the 11th of Febuary, 2020 (Tuesday) from 6am to 10am. We then compare the
total outbound and inbound commuters output by the algorithm with the census
data. The results are displayed in Figure 10, including only those census
regions which clearly correspond to one or more SA-2 regions. Moreover, cell
coverage is not necessarily representative of an individual’s exact location
at any given moment, and data between neighbouring regions may be somewhat
mixed.111111The code used to generate the comparisons shown here, along with
the corresponding initial conditions, is in the folder ‘southland-comparison’.

Figure 10: Stats NZ Commuter View vs. movements estimated with the “exact” and
“approximate” algorithm for Southland commuters. We used data from 11th
February, for the five hours from 6am to 10am. SA-2 regions which do not
correspond to a single commuter region are discarded in this analysis.

While the algorithm appears to capture the essential features of the inbound
commuting trends, with gathering concentrated within the main metropolis, the
outbound inference fares significantly less well, with outbound commuters
significantly underrepresented when compared to regions within the census data
with high traffic. This comparison is strictly qualitative and “one off”, and
self-declared travel within inner-city regions may correspond to a change in
cell-coverage regions, and vice-versa.

We also note that the ‘Commuter View’ tool has since been re-branded to
‘Commuter Waka’, and now also incorporates data from the 2018 census [10].
However, due to a particularly low response rate of 81.6% to the 2018 census
[11], we choose to test our algorithm against the older data set - based upon
the responses to the 2013 census only - which had a much higher response rate.
It is hoped that better quality data will become available in future for more
thorough verification testing.

##  9 Summary

This Working Paper describes and re-implements two possible approaches to
using count data to impute meso-scale human movement patterns. We investigated
the validity of the likelihood assumed and improved the minimization method.
At this point it can analyse data for large fractions of the country (e.g.
$\sim$800 out of $\sim$2000 SA-2 regions in a single run) via a
computationally efficient and publicly available code.

The algorithm demonstrates qualitative agreement with simulated and real-world
data. The actual numerical counts of people moving from one region to another
come with some uncertainty, stemming from the fact that the problem is highly
degenerate. In particular, we occasionally find estimated values of
$\mathbf{s}$, $\bm{\pi}$ $\mathbf{M}$ and $\beta$ that have a higher
likelihood than the “true” solution, but nevertheless differ from it.
Moreover, the model used here computes large numbers of “point to point”
journeys, but in real world scenarios residents of two well-separated blocks
may be more likely to interact in some third block, to which they both travel
during the course of day.

That said, we can see a number of ways in which these approaches could be
improved, and our implementations are publicly available. The algorithms could
prove useful when averaged data is extracted from the outputs, such as
estimates of mean distances travelled per day. Such averaged quantities may be
more reliable than estimates of individual point-to-point journeys. The codes
are therefore probably best regarded as providing order-of-magnitude estimates
which may be of use in sensitivity testing complex infection propagation
models, but should not be seen as yielding precise quantitative predictions.

While improvement is still needed, our work may have important applications in
many areas relating to disease outbreak and containment. Namely,
identification of the areas with highest gathering statistics could help to
inform the most effective locations for lockdown boundaries, while a better
understanding of common transit routes could help to identify high risk sub-
regions outside of the most densely populated commercial and residential hubs.
Finally, outputs from this algorithm may serve as useful inputs to more
complex models of disease propagation.

Specific directions for future work might include:

  * •

Adding a more nuanced distance metric, including driving distance, rather the
centroid to centroid Euclidean distance.

  * •

Considering a more complex penalty function, e.g. $\exp(-\beta d-\alpha
d^{2})$.

  * •

Improving the quality of the data set. In particular, a count of cell phones
in a block that were present in the previous hour would allow separate
estimations of the $s_{i}$ and $\pi_{i}$ and would fix the diagonal elements
of $\mathbf{M}$, but would likely raise few privacy concerns.

  * •

Improving validation testing against census data, or traffic flow information
for urban regions.

  * •

Fitting $\bm{\pi}$, $\mathbf{s}$, and (possibly) $\beta$ to census data rather
than count data. One would have to justify the continued use of these values
during periods of modified behaviour, such as when travel restrictions are in
place.

  * •

Developing an improved travel simulator to better test the model against a
full realisation of movement patterns in areas of interest.

  * •

Properly accounting for the fact that in most realistic datasets the
transition probabilities will be time dependent, varying over the course of
the (working) day.

Finally, we emphasise that this overall problem is essentially an imputation
exercise. Any results obtained with them are estimates, and any model that
uses them as inputs should be interpreted accordingly.

## Biographical Note

This work was undertaken in response to the Covid-19 emergency and was largely
completed while New Zealand was at lockdown levels 3 and 4. At the time the
work was undertaken the authors were all members of the theoretical cosmology
group at the University of Auckland.

## References

  * [1] Yasunori Akagi, Takuya Nishimura, Takeshi Kurashima, and Hiroyuki Toda.  A Fast and Accurate Method for Estimating People Flow from Spatiotemporal Population Data.  In Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence, pages 3293–3300, Stockholm, Sweden, July 2018. International Joint Conferences on Artificial Intelligence Organization.
  * [2] Richard H. Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu.  A Limited Memory Algorithm for Bound Constrained Optimization.  SIAM Journal on Scientific Computing, 16(5):1190–1208, September 1995.  Publisher: Society for Industrial and Applied Mathematics.
  * [3] A. P. Dempster, N. M. Laird, and D. B. Rubin.  Maximum likelihood from incomplete data via the em algorithm.  Journal of the Royal Statistical Society. Series B (Methodological), 39(1):1–38, 1977.
  * [4] Chuong B. Do and Serafim Batzoglou.  What is the expectation maximization algorithm?  Nature Biotechnology, 26(8):897–899, Aug 2008.
  * [5] Nir Friedman.  The Bayesian Structural EM Algorithm.  arXiv e-prints, page arXiv:1301.7373, January 2013.
  * [6] Siu Kwan Lam, Antoine Pitrou, and Stanley Seibert.  Numba: A llvm-based python jit compiler.  In Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, LLVM ’15, New York, NY, USA, 2015. Association for Computing Machinery.
  * [7] Xiao-Li Meng and Donald B. Rubin.  Maximum likelihood estimation via the ecm algorithm: A general framework.  Biometrika, 80(2):267–278, 1993.
  * [8] Stats NZ.  Census: Household composition by tenure of household and access to telecommunication systems, for households in occupied private dwellings 2018.  http://nzdotstat.stats.govt.nz/WBOS/Index.aspx?DataSetCode=TABLECODE8462.  Accessed: 2023-01-18.
  * [9] Stats NZ.  Commuter view.  https://www.stats.govt.nz/tools/commuter-view.  Accessed: 2020-07-13.
  * [10] Stats NZ.  Commuter waka.  https://commuter.waka.app.  Accessed: 2023-01-19.
  * [11] Stats NZ.  Post-enumeration survey: 2018.  https://www.stats.govt.nz/information-releases/post-enumeration-survey-2018/.  Accessed: 2023-01-19.
  * [12] Stats NZ.  Statistical area 2 2018 v1.0.0.  https://datafinder.stats.govt.nz/layer/92212-statistical-area-2-2018-generalised/.  Accessed: 2020-06-05.
  * [13] Stéfan van der Walt, S. Chris Colbert, and Gaël Varoquaux.  The NumPy Array: A Structure for Efficient Numerical Computation.  Comput. Sci. Eng., 13(2):22–30, 2011.
  * [14] Pauli Virtanen et al.  SciPy 1.0–Fundamental Algorithms for Scientific Computing in Python.  Nature Meth., 2020.


11institutetext: Rapid7, Boston MA 02114, USA
22institutetext: Drexel University, Philadelphia PA 19104, USA
22email<EMAIL_ADDRESS>
# Winning the Ransomware Lottery ††thanks: Funded in part by the Auerbach
Berger Chair in Cybersecurity held by Spiros Mancoridis, at Drexel University

A Game-Theoretic Approach to Preventing Ransomware Attacks

Erick Galinkin 11 2 2 0000-0003-1268-9258

###### Abstract

Ransomware is a growing threat to individuals and enterprises alike,
constituting a major factor in cyber insurance and in the security planning of
every organization. Although the game theoretic lens often frames the game as
a competition between equals – a profit maximizing attacker and a loss
minimizing defender – the reality of many situations is that ransomware
organizations are not playing a non-cooperative game, they are playing a
lottery. The wanton behavior of attackers creates a situation where many
victims are hit more than once by ransomware operators, sometimes even by the
same group. If defenders wish to combat malware, they must then seek to remove
the incentives of it. In this work, we construct an expected value model based
on data from actual ransomware attacks and identify three variables: the value
of payments, the cost of an attack, and the probability of payment. Using this
model, we consider the potential to manipulate these variables to reduce the
profit motive associated with ransomware attack. Based on the model, we
present mitigations to encourage an environment that is hostile to ransomware
operators. In particular, we find that off-site backups and government
incentives for their adoption are the most fruitful avenue for combating
ransomware.

###### Keywords:

Security Malware Economics Ransomware Incentives Backups

##  1 Introduction

Ransomware is a family of malware that encrypts files on a system and demands
payment for the ability to decrypt these files. Although proof of concept
ransomware has existed since at least 1996 [35], modern ransomware tactics
result from CryptoLocker’s revolutionary use of Bitcoin for payment [14]. This
innovation has allowed ransomware actors to perpetrate increasingly
sophisticated attacks, including the 2017 WannaCry attack [16] – an attack
whose effects, according to ransomware payment tracker
Ransomwhere111https://ransomwhe.re are still being felt today. We have seen a
pivot in targeting, from the wanton use of exploit kits and watering hole
attacks that largely affected end users to the current increase in enterprise
victims [22] by way of malicious loaders and initial access brokers [12].

The threat of ransomware grows larger year after year, with a spate of recent
attacks including on the Colonial pipeline [18] and the Kaseya supply chain
attack [23] demonstrating the devastation and real-world impact of the issues.
The Ransomware Task Force report [9] identifies the goal of disrupting the
ransomware business model as an important goal. This goal is uniquely
important, since ransomware is so often an attack of opportunity – akin to a
mugging or kidnapping – and not the sort of highly-targeted attack that is
often expected from sophisticated adversaries. We frame the problem in a new
way, as the attacker is not playing a single game against a single defender.
Rather, attackers seek to find vulnerable victims wherever they may be, and so
instead of playing a game with attackers, we view the problem from the
attacker point of view. To this end, we suggest that defenders should consider
the problem of ransomware and ransomware payments in particular as analogous
to an attacker playing a lottery instead of a strategic game between equals.

##  2 Related Work

In recent years, considerable research has been done on the game theory of
ransomware payments. The earliest relevant work on the topic appears to be by
Spyridopoulos et al. [27], who found a Nash equilibrium balancing potential
costs of mitigation with the cost of a successful attack. Leveraging
epidemiologically-inspired models of malware spread, this work considered the
equilibria of available defender strategies. The game is constructed under a
unified proliferation model, with infection, immunization, and disinfection
rates that informed the strategies of the players. These player’s payoffs were
then computed for a set of strategies given the parameters controlled by the
attacker and the defender – the infection rate, patch rate, removal rate, and
the rate of both patching and removal. Spryidopoulos et al.’s work informed
defenders how to approach ransomware worm attacks and defined the optimal
strategy for the defender.

The work of Laszka et al. [13] was the first to consider the economics of
ransomware using models that reflect the similarity of ransomware to
kidnapping and ransom. They developed an economic model of the interaction
between attackers and victim organizations, and studied that model to minimize
the economic impact to those organizations. Primarily, the work focused on the
cost-benefit of investing in backup solutions, a recommendation that is still
widely regarded as the best way to prepare for ransomware attacks [9]. Laszka
et al. also showed how coordinated backup investments can deter ransomware
attackers in particular – a novel insight in the literature. Our work borrows
from their recommendations and builds on this existing literature, but we
differ in our approach to the game-theoretic model.

Caporusso et al. [5] also built upon the kidnap and ransom literature,
leveraging a negotiation model represented as an extensive-form game. This
work dealt with ransomware in cases where renegotiation of the ransom is
possible, a surprisingly common phenomenon that has been seen with some
ransomware operators [17] – though other ransomware operators refuse to
negotiate. Caporusso et al. identified the post-attack dynamics between the
human victim and the human ransomware operator, acknowledging that there are
substantial human factors outside of ransom negotiation to be made in the
decision making process.

Cartwright et al. [6] grappled with the question of whether or not to pay a
ransom at all. Their work largely built upon the earlier paper of Laszka et
al. and framed the problem of ransomware under the lens of kidnap and ransom.
It did so by building upon two existing kidnapping models, those of Selten
[25], and Lapan and Sandler [11]. The Selten model informed the optimal ransom
to be set by the attacker, while the model of Lapan and Sandler aided in
deciding whether or not victims should take action to deter the kidnapping in
the first place. In contrast to this work, we present a novel approach to the
game and develop a model under a differing set of assumptions.

##  3 Probability and Lotteries

In common parlance, “lottery” typically refers to a form of gambling where a
player purchases a ticket at some nominal cost with a fixed set of different
numbers. Then, another set of numbers with the same size is drawn at random
without replacement. After this draw, some reward that confers some amount of
utility may be given depending on how many numbers in the randomly drawn set
match the set on the purchased ticket.

Mathematically, we can formalize a lottery as follows: Let $X$ be a set of
prizes, $X=\\{x_{1},...,x_{n}\\}$, that confers some utility. From this set of
prizes, we define a lottery $L=\\{p_{i},...,p_{n}\\}$ over the set of prizes
such that for each $x_{i}\in X$, there is a corresponding $p_{i}\geq 0$, and
$\sum_{i=1}^{n}p_{i}=1$. There is also some cost $c\geq 0$ to enter the
lottery. Then, for each of the prizes, there is some utility $u(x_{i})$ that
the agent derives from receiving that prize, and their expected utility over
the lottery is then $\sum_{i=1}^{n}p_{i}u(x_{i})-c$. In the ransomware
context, a prize $x$ corresponds to a payment to a ransomware operator, and
$p$ is the probability that a victim will pay that amount.

The optimal ransom value for $x$ has been explored in other work [6] so we
instead deal with the binary probability that a victim will pay or not pay,
assuming that the optimal ransom value is set. In our ransomware lottery, we
thus define 2 probabilities: $p_{\text{win}}$, when a victim pays a ransom and
$p_{\text{lose}}=1-p_{\text{win}}$, when a victim does not. For simplicity in
this initial model, we incorporate the probability that the attack is not
successful into $p_{\text{lose}}$. There is, as mentioned, also some small
cost $c$ associated with launching the ransomware attack.

Conveniently for ransomware operators, $c$ is quite small, and
$x_{\text{win}}$ can be quite large, as we discuss in Section 4. By contrast,
$x_{\text{lose}}=0$, since there is no chance that ransomware operators will
have to pay more than the cost to launch the attack – the victim will simply
ignore the attack because they do not value the information which has been
ransomed or have some mitigation such as those outlined in Section 5. In
total, this means that the game played, from the perspective of ransomware
operators, is as follows:

$\displaystyle L$ $\displaystyle=\\{p_{\text{win}},p_{\text{lose}}\\}$
$\displaystyle X$ $\displaystyle=\\{x_{\text{win}},0\\}$

and therefore, the expected utility for a single successful attack is:

$\displaystyle E[u(x)]$
$\displaystyle=\sum_{i=\\{\text{win},\text{lose}\\}}p_{i}(x_{i}-c)$
$\displaystyle=(p_{\text{win}}(x_{\text{win}}-c))+(p_{\text{lose}}(0-c))$
$\displaystyle=p_{\text{win}}x_{\text{win}}-(p_{\text{win}}c+p_{\text{lose}}c)$
$\displaystyle=p_{\text{win}}x_{\text{win}}-c$ (1)

Since $x_{\text{lose}}=0$ and $p_{\text{lose}}=1-p_{\text{win}}$, for the sake
of simplicity and readability, we use $x$ and $p$ in the remainder of the
paper to represent the case when a victim pays. We can see from Equation 1
that ransomware operators are incentivized to continue operating for as long
as the value of $px>c$, since they will profit from each attack, on average.
Research by Kaspersky Labs [10] shows that 56% of ransomware victims pay the
ransom to restore access to their data. At this rate of payment, the cost of
an average ransomware attack would need to be 1.7857 times – nearly double –
the optimal payment to remove the incentive.

We can see that probabilistically, this is equivalent to betting on a biased
coin flip. Since $E[u(x)]$ is a function of the random variable $x$, it is
itself a random variable, which we denote $Y$. Given a cost to make a bet $c$,
we flip a biased coin with win probability $p$ and receive payout $x$ at that
rate. Let $b$ be the amount of capital available to the bettor – our attacker
– and let $b>c$. We initialize $b_{0}$ to be the amount of capital available
before any bets are cast and $b_{i}$ the available capital to the bettor at
trial $i$. Then after the first trial, our possible values for $b_{1}$ are
$b_{1}=b_{0}-c$ or $b_{1}=b_{0}-c+x$. Our expected value of
$b_{1}=(b_{0}-c)+px$, as in Equation 1.

By the linearity of expectation, our expected bank at trial $k$ is:

$b_{k}=b_{0}+E[Y_{k}]=b_{0}+k(px-c)$

We can see that if $px>c$, then the expected value of each trial is positive,
and so for the player making the bet,

$\lim_{k\rightarrow\infty}E[Y_{k}]=k(px-c)=\infty$ (2)

This suggests that any player who can participate in the game is highly
incentivized to play as many rounds as possible, since the potential payoff is
infinite. Note that this expected value only holds in an idealized world with
infinite money and no law enforcement, so it does not capture the intricate
relationships of the real world. It does, however, demonstrate that since the
expectation is not finite, there is no optimal stopping time. Therefore, there
is no incentive for any attacker to ever stop conducting ransomware attacks
when $px-c$ is reasonably large.

To demonstrate this, we construct three simple simulations, shown in Figure 1.
We set our payout value $x=170404$ and cost $c=4200$ based on analysis in
Section 4. Then, for three different values of $p$: 0.1, 0.3024, and 0.5, we
run 1000 trials. With probability $p$, the player receives value $x-c$, and
with probabiltiy $1-p$, the player receives value $-c$. We can see that
overall, the accumulated value is linear with respect to $p$, as we would
expect from Equation 1.

Figure 1: Plot of simulation demonstrating accumulated utility at p=0.1,
p=0.3024, and p=0.5

##  4 Paying to Play

The cost of running a ransomware attack is very opaque and highly variable.
Some cybercriminal organizations are sophisticated operations that develop
their malware in-house [31]. These organizations have software development
lifecycles, version control, testing, and pay staff to perform all of these
functions. Other organizations simply purchase ransomware-as-a-service [15]
(RaaS) or piece together their arsenal from so-called darknet markets. A 2017
study [30] found that prices ranged from $0.50 to $3,000 for ransomware
products, at a median price of $10.50. In contrast to these prices, most RaaS
providers take a percentage of the ransom, rather than providing an executable
for a flat fee.

In order to infect an endpoint with ransomware, however, one needs to gain
initial access. Furthermore, most ransomware operators leverage a loader – a
small program designed to install another malware on a target system – to
actually get the ransomware onto the endpoint. Nearly all ransomware variants
[20] rely on phishing, commodity malware, exploit kits, and vulnerable
services – particularly the remote desktop protocol – to deliver their
malware. This factors in to the overall cost of operation, but is challenging
to estimate, since cybercriminals are not forthcoming with this information. A
technical report issued by Deloitte [1] found the cost of initial access to be
between $70 and $400 per 1000 machines depending on geographic region, and the
cost of a loader to range from $3 to $4,000, depending on functionality. The
United States demanded the highest fee for an initial access at $400. At this
time, the US is also the nation which demands the highest ransoms, and so in
the interest of creating a conservative but accurate estimate, we use this
number. The highest average monthly cost of a loader was $800, which is the
figure we use moving forward. We thus estimate the cost of an attack at
$c=3000+400+800=4200$.

This cost of $4,200 means at at a payment rate of $p=0.56$, the minimal ransom
to turn a profit is $7,500. However, this payment rate is too large, since it
assumes that the attack has been successful. According to Sophos [26], only
54% of attacks actually encrypt data. Given that a successful attack is a
precondition for being a paying victim, the joint probability of the attack
being successful and the ransom being paid, which we defined in Equation 1 as
$p_{\text{win}}$ is the product of these two probabilities. Our joint
probability for a successful attack where the victim pays the ransom is
therefore:

$p=P(\text{paid}|\text{success})\cdot P(\text{success})=0.56\cdot 0.54=0.3024$

This suggests that at a cost of $4,200, per attack the minimal ransom an
attacker must request to remain profitable is $13,888.89. As of March 2021,
the average value of ransomware a payout for a compromised organization was
$312,493 [8], around 22 times the minimal value needed to incentivize the
attacks. We note that other estimates, such as those by Sophos [26] are a more
modest $170,404 for mid-sized organizations in the United states, a value
which is still around 12 times the minimum to create positive expected value
for these attacks. We treat these as a “reasonable average range” in our
subsequent analysis.

There are three variables in this problem that may disincentivize the
perpetration of ransomware attacks:

  1. 1.

Lowering the value of the payments

  2. 2.

Increasing the cost of operating ransomware

  3. 3.

Decreasing the probability of payment

We discuss the feasibility of using each of these three variables to
disincentivize ransomware attacks in turn.

###  4.1 Lowering the Value of Payments

Today, there are few options for lowering the value of a payment. Since nearly
all payments for ransomware are rendered in cryptocurrency, a steep decline in
the value of cryptocurrency or the inability to exchange it for other goods or
services would remove the effective value of a successful attack. To date,
some proposals have been made to ban [7], or regulate cryptocurrencies [19,
24], though the effect of these bans and proposed regulations on the price of
cryptocurrency remains to be seen. Moreover, even if cryptocurrency were
regulated into obsolescence, ransoms could be paid in gift cards or other hard
to track currency equivalents. This suggests that lowering the value of
payments is not a viable path for removing the incentive.

###  4.2 Increasing Costs

The onus for increasing costs falls on the ransomware developers and operators
themselves, and so there is likely a cost ceiling. If the marketplace
efficiencies of initial access brokers and ransomware-as-a-service were
removed entirely, the cost of conducting an attack would be the cost of
development plus the cost of deployment and maintenance of the infrastructure.
This would require more technical skill and initial investment than relatively
low-skill ransomware operators would be capable of, but after the initial
investment, would likely cost less per-attack than the $3,000 high-end figure
from [30]. This may, on balance, reduce the overall prevalence of malware
attacks. However, this would also require the takedown of nearly all darknet
marketplaces. Despite a number of high-profile takedowns, ransomware continues
to flourish on these marketplaces. Thus, the options for increasing costs to
operators are also limited.

###  4.3 Decreasing Payment Probability

Since the probability of payment is the one thing out of the control of the
attackers, it stands to reason that it is where defenders can exercise the
most control. In our model, decreasing the probability of a successful attack
that gets paid linearly reduces the expected value of an attack. This means
that organizations have two options available to them to reduce an attack’s
expected value. Decreasing the success of launched attacks will prevent the
victim having to decide whether or not to pay the ransom in the first place.
Assuming an attack is successful, decreasing the chance that the ransom is
paid will also reduce the attacker’s value.

Given our average payout value range of $x=[170,404,312493]$, the expected
value of an attack at current payment rates is in the range
$[47,300.17,170,798.08]$. A 50% reduction in probability of payout to $p=0.28$
against a cost of $c=4200$, with attack success rates held equal yields an
expected value range of $[21565.08,43048.94]$ – an amount that a would-be
ransomware operator could make as a software engineer in Europe [21] instead
of perpetrating ransomware attacks. Given the financial motivation of most
ransomware operators [2], it stands to reason that a comparable salary is a
perfectly substitutable good for rational actors. To eliminate profit
entirely, assuming current attack success rates and sufficient economies of
scale, payment probability would need to decrease to 2.489% on the high-end of
average payments and 4.564% on the low-end of payments – a dramatic reduction
from today’s payment rates.

Despite that “break-even” probability, ransomware operators are likely to turn
to some other income stream before profits hit zero due to law enforcement
activities surrounding cybercrime. In particular, the US Federal Bureau of
Investigations and the UK National Cyber Security Centre have pursued
cybercriminals abroad [28], indicting and sanctioning ransomware operators.
However, in order to drastically reduce the payout rate of ransomware,
organizations will need to have a reason not to pay the ransoms.

##  5 Lowering the Stakes

In order to lower the probability of payment and create an environment where
attackers are not incentivized to continue launching ransomware attacks,
victims must be incentivized not to pay the ransom. An effective strategy for
lowering the probability of payment ultimately consists of one where the
victim’s options for restoration are meaningfully less costly than paying the
ransom. Considerable work has been done on quantifying these differences and
we point to the article by Cluley [8] for details, as the specific rates will
differ from organization to organization. Since the use of ransomware is
illegal, there are external, non-financial mechanisms for reducing attacker
incentives such as arrest, seizure of assets, indictment, and sanctions. We do
not address these mechanisms in our framework and reserve their impact for
future work.

In order to reduce attacker incentives, we consider the potential impact of
four commonly discussed strategies:

  1. 1.

Decreasing Attack Success

  2. 2.

Cyber Insurance

  3. 3.

Use of Decrypters

  4. 4.

Off-Site Backups

###  5.1 Decreasing Attack Success

Decreasing attack success is the goal of any organizational information
security program. The success of attacks has myriad factors, ranging from
human factors such as insider threats and phishing to software vulnerabilities
and misconfigurations. Modern antivirus technologies can assist in catching
the loaders that often deliver the ransomware, and some endpoint security
solutions can even detect exploitation of vulnerabilities. In addition,
training programs for phishing emails and advising customers not to open
attachments from unknown senders are widely used to attempt to mitigate these
attacks. A comprehensive listing of ways to reduce an organization’s attack
surface is out of the scope of this paper, but a 2020 report by Deloitte and
the Financial Services Information Sharing and Analysis Center [4] showed that
on average, 10% of an organization’s information technology budget –
approximately 0.2% of company revenue – is dedicated to cybersecurity. In
light of the increasing threats associated with ransomware, this amount may
not be sufficient to reduce the probability that an attack is successful.

The figure in Equation 1 only holds for cases where a ransomware infection has
been successful and does not account for failed attacks – only payments.
Reducing the incidence of these attacks through other means such as the use of
application allowlists, strong spam filters, protection of exposed ports and
services, and other well-known security hygiene methods can serve to reduce
the success of these attacks. Since the cost to an attacker is undertaken
whether or not the attack is successful, the failure of these attacks will
discourage these attackers. In order to isolate the influence of payment
probability, our analysis assumed that all attacks are successful – a naive
assumption that suggests the 1.5% payout probability derived in Section 4.3 is
the probability of payment overall, not merely the conditional probability of
payment given a successful attack.

###  5.2 Cyber Insurance

Cyber insurance is a strategy that is often mentioned as an organizational
solution in the context of ransomware. This can help to protect businesses
from the cost of ransomware attacks, covering the cost to restore encrypted
data. However, in cases where cyber insurance alleviates the burden to
victims, attackers are still paid, doing nothing to remove the incentives
surrounding ransomware. Consequently, from an attacker incentive perspective,
cyber insurance does nothing to alleviate the overall problem of ransomware.

###  5.3 Use of Decrypters

The use of decrypters is a significant way to allow victims to ignore the
effects of ransomware. Although decrypters for some of the most popular
strains of ransomware today are not available, organizations like No More
Ransom!222https://www.nomoreransom.org offers free decrypters for more than
150 families of ransomware. Widespread knowledge of these utilities and
increased investment by security researchers on developing these utilities
could allow victims to decrypt their own files without paying a ransom. Note
that when decrypters become available or kill-switches as seen in WannaCry
[16] shut down operations, ransomware operators will patch their malware [3]
to continue operations.

###  5.4 Off-Site Backups

The most commonly proposed solution for organizations to avoid the impacts of
ransomware and confidently be able to not pay a ransom is the use of off-site
backups. An off-site backup can be used to restore systems to pre-ransomware
configurations and tends to cost significantly less than paying the ransom.
Research by Wood et al. [34] acknowledges the difficulties of backup
deployments. Although they develop their recovery from a disaster preparedness
perspective, their cost estimates show that both cloud-based and colocation
for backups can allow for high uptime at a fraction of the cost associated
with paying a ransom. Additionally, having a backup that allows for
restoration reduces the cost to remediate possible residual traces of the
attacker, reduces time to remediate, and mitigates much of the reputational
damage associated with paying a ransom.

###  5.5 Impact of Mitigations

The aforementioned approaches may allow victims to choose not to pay, but as
Cartwright et al. [6] demonstrate, victims will have different willingness to
pay given some set ransom. This willingness to pay depends on the size of the
ransom and therefore encourages the victim to mitigate the attack. When
victims pay, they usually – though not always [26] – get their files back, a
factor which discourages paying. However, there is some cost to deterrence,
and if that is too high, the victim will instead accept their chances of being
infected.

There are also factors at play external to the relationship between the cost
of a ransom versus the cost of mitigation. For example, in the United States,
ransom payments can be written off [33] as “ordinary, necessary, and
reasonable” expenses for tax purposes. This factor actually incentivizes
victims to pay, and discourages additional investments into mitigation.
Wheeler and Martin [32] point out that in the current regulatory environment
of the United States, there is a misalignment between public interests to
discourage ransomware and private interests to recover data and resume
operations at the lowest cost. We conclude then, that government and
regulatory organizations interested in preventing ransomware should create
financial incentives for organizations and individuals to invest in backups
that allow for ransoms not to be paid. Further, policy solutions to change the
tax incentives associated with paying ransoms could be pursued to improve the
chance that companies will invest in security technologies.

##  6 Conclusion

Ransomware remains a significant problem in the world, and our analysis
demonstrates why – there is effectively unlimited incentive to use ransomware.
Since the cost is relatively low and the potential payouts are high,
financially-motivated actors are encouraged to pursue this line of attack.
Additionally, the victims of successful attacks are more likely to pay than
not for a variety of factors, including the ability to write-off the ransom as
a business expense.

If we wish to eliminate the threat of ransomware, we cannot attack the market
itself, as the actors are aware that their actions are illegal but have
accepted that risk. Instead, we must see that attackers are engaged in a
simple game where they do not need to account for the strategies of their
victims. Where defenders have power to affect ransomware is largely on the
front of actually paying the ransoms.

We outlined a handful of commonly-discussed solutions and conclude that off-
site backups remain the most effective way to ignore the impact of ransomware
attacks. In order to encourage organizations to pursue these policies, we
conclude that governmental and regulatory organizations will need to provide
incentives for organizations to invest in these backup solutions. Short of
encouraging these solutions and allowing victims not to pay ransoms, we can
reasonably expect the ransomware threat to continue to grow.

The model used here leveraged a probabilistic model and expected utility
theory to identify incentives and explore the security impacts of those
incentives. In future work, we seek to explore a more realistic model of the
risk behaviors these attackers and defenders exhibit based on their subjective
beliefs. Furthermore, there are meaningful non-financial mechanisms such as
those mentioned in Section 5, and inclusion of those mechanisms would require
a more complex model. This could be done by representing uncertainty via
cumulative prospect theory [29], as has been done in the economic literature.
In particular, there is a significant amount of uncertainty on the part of
attackers about whether or not an attack will be successful. Similarly, there
is significant uncertainty for defenders about how, when, and where they will
be attacked. By representing the choice under uncertainty more richly than in
an expected utility model, we may better model the true behaviors of attackers
and defenders.

## References

  * [1] Analytics, D.T.I..: Black-market ecosystem: Estimating the cost of “pwnership”. Deloitte Technical Report (2018), https://www2.deloitte.com/us/en/pages/risk/articles/vigilant-threat-studies-deloitte-us.html
  * [2] Anderson, R.: Security engineering: a guide to building dependable distributed systems. John Wiley & Sons (2020)
  * [3] Arghire, I.: “patched” wannacry ransomware has no kill-switch. SecurityWeek (2017), https://www.securityweek.com/patched-wannacry-ransomware-has-no-kill-switch
  * [4] Bernard, J., Nicholson, M.: Reshaping the cybersecurity landscape. Deloitte Technical Report (2020), https://www2.deloitte.com/us/en/insights/industry/financial-services/cybersecurity-maturity-financial-institutions-cyber-risk.html
  * [5] Caporusso, N., Chea, S., Abukhaled, R.: A game-theoretical model of ransomware. In: International Conference on Applied Human Factors and Ergonomics. pp. 69–78. Springer (2018)
  * [6] Cartwright, E., Hernandez Castro, J., Cartwright, A.: To pay or not: game theoretic models of ransomware. Journal of Cybersecurity 5(1), tyz009 (2019)
  * [7] Clark, M.: What we know about china’s cryptocurrency crackdown. Vox (2021), https://www.theverge.com/2021/6/23/22544367/china-crypto-crackdown-bitcoin-mining-sichuan-ban-hydro-cryptocurrency-trading
  * [8] Cluley, G.: Average ransomware payouts shoot up 171% to over $300,000. Tripwire – The State of Security (2021), https://www.tripwire.com/state-of-security/featured/average-ransomware-payouts-shoot-up/
  * [9] Force, R.T.: Combating ransomware (2021)
  * [10] labs, K.: Consumer appetite versus action: the state of data privacy amid growing digital dependency. Kaspersky Consumer IT Security Risks Report 2021 (2021), https://media.kasperskydaily.com/wp-content/uploads/sites/92/2021/03/16090300/consumer-appetite-versus-action-report.pdf
  * [11] Lapan, H.E., Sandler, T.: To bargain or not to bargain: That is the question. The American Economic Review 78(2), 16–21 (1988)
  * [12] Larson, S., Blackford, D., G, G.: The first step: Initial access leads to ransomware. Proofpoint Threat Insight (2021), https://www.proofpoint.com/us/blog/threat-insight/first-step-initial-access-leads-ransomware
  * [13] Laszka, A., Farhang, S., Grossklags, J.: On the economics of ransomware. In: International Conference on Decision and Game Theory for Security. pp. 397–417. Springer (2017)
  * [14] Liao, K., Zhao, Z., Doupé, A., Ahn, G.J.: Behind closed doors: measurement and analysis of cryptolocker ransoms in bitcoin. In: 2016 APWG symposium on electronic crime research (eCrime). pp. 1–13. IEEE (2016)
  * [15] Meland, P.H., Bayoumy, Y.F.F., Sindre, G.: The ransomware-as-a-service economy within the darknet. Computers & Security 92, 101762 (2020). https://doi.org/https://doi.org/10.1016/j.cose.2020.101762, https://www.sciencedirect.com/science/article/pii/S0167404820300468
  * [16] Mohurle, S., Patil, M.: A brief study of wannacry threat: Ransomware attack 2017\. International Journal of Advanced Research in Computer Science 8(5), 1938–1940 (2017)
  * [17] Monroe, R.: How to negotiate with ransomware hackers. The New Yorker (2021), https://www.newyorker.com/magazine/2021/06/07/how-to-negotiate-with-ransomware-hackers
  * [18] Morrison, S.: How a major oil pipeline got held for ransom. Vox (2021), https://www.vox.com/recode/22428774/ransomeware-pipeline-colonial-darkside-gas-prices
  * [19] Nabilou, H.: How to regulate bitcoin? decentralized regulation for a decentralized cryptocurrency. International Journal of Law and Information Technology 27(3), 266–291 (2019)
  * [20] Networks, P.A.: Ransomware threat report, 2021. Palo Alto Networks Technical Report (2021), https://www.paloaltonetworks.com/resources/research/unit42-ransomware-threat-report-2021
  * [21] Orosz, G.: The trimodal nature of software engineering salaries in the netherlands and europe. Pragmatic Engineer (2021), https://blog.pragmaticengineer.com/software-engineering-salaries-in-the-netherlands-and-europe/
  * [22] O’Gorman, B., Wueest, C., O’Brien, D., Cleary, G.: Symantec internet security threat report. Symantec Corp., Mountain View, CA, USA, Tech. Rep (2019)
  * [23] Press, A.: Scale, details of massive kaseya ransomware attack emerge. NPR (2021), https://www.npr.org/2021/07/05/1013117515/scale-details-of-massive-kaseya-ransomware-attack-emerge
  * [24] Schaupp, L.C., Festa, M.: Cryptocurrency adoption and the road to regulation. In: Proceedings of the 19th Annual International Conference on Digital Government Research: Governance in the Data Age. pp. 1–9 (2018)
  * [25] Selten, R.: Models of strategic rationality, vol. 2. Springer Science & Business Media (2013)
  * [26] Sophos: Sophos state of ransomware 2021. Sophos Technical Report (2021), https://secure2.sophos.com/en-us/medialibrary/pdfs/whitepaper/sophos-state-of-ransomware-2021-wp.pdf
  * [27] Spyridopoulos, T., Maraslis, K., Mylonas, A., Tryfonas, T., Oikonomou, G.: A game theoretical method for cost-benefit analysis of malware dissemination prevention. Information Security Journal: A Global Perspective 24(4-6), 164–176 (2015)
  * [28] Tidy, J.: The ransomware surge ruining lives. BBC (2021), https://www.bbc.com/news/technology-56933733
  * [29] Tversky, A., Kahneman, D.: Advances in prospect theory: Cumulative representation of uncertainty. Journal of Risk and uncertainty 5(4), 297–323 (1992)
  * [30] Unit, C.B.T.A.: Dark web ransomware economy growing at an annual rate of 2,500%. Carbon Black Threat Research (2017), https://www.carbonblack.com/2017/10/11/dark-web-ransomware-economy-growing-annual-rate-2500/
  * [31] U.S. Attorney’s Office, Western District of Washington: High-level organizer of notorious hacking group fin7 sentenced to ten years in prison for scheme that compromised tens of millions of debit and credit cards (2021)
  * [32] Wheeler, T., Martin, C.: Should ransomware payments be banned? The Brookings Institute Tech Stream (2021), https://www.brookings.edu/techstream/should-ransomware-payments-be-banned/
  * [33] Wood, R.: Garmin hack’s $10m ransom payment, $10m tax deduction. Forbes (2020), https://www.forbes.com/sites/robertwood/2020/07/27/garmin-hacks-10m-ransom-payment-10m-tax-deduction/?sh=4452ae4712c5
  * [34] Wood, T., Cecchet, E., Ramakrishnan, K.K., Shenoy, P.J., van der Merwe, J.E., Venkataramani, A.: Disaster recovery as a cloud service: economic benefits & deployment challenges. HotCloud 10, 8–15 (2010)
  * [35] Young, A., Yung, M.: Cryptovirology: extortion-based security threats and countermeasures. In: Proceedings 1996 IEEE Symposium on Security and Privacy. pp. 129–140. IEEE (1996)