\documentclass{article}
\usepackage[utf8x]{inputenc}
\usepackage{amsmath, amsthm}
\usepackage{algorithm}% http://ctan.org/pkg/algorithms
\usepackage{algpseudocode}% http://ctan.org/pkg/algorithmicx
\usepackage{graphicx}
\usepackage{microtype}
\usepackage{verbatim}

\DeclareMathOperator*{\argmax}{argmax}
\newcommand{\eqdef}{\mathbin{\stackrel{\rm def}{=}}}
\providecommand{\e}[1]{\ensuremath{\times 10^{#1}}}
\newcommand{\given}{\,|\,}
\newtheorem{theorem}{Theorem}
\newtheorem{proposition}[theorem]{Proposition}

\title{Hawkes Contagion model}
\author{Ben Green \and Thibaut Horel \and Andrew Papachristos}
\date{September 2015}

\begin{document}

\maketitle

\section{Contagion Model}

We model the contagion of violence using a multi-dimensional Hawkes process
over the co-offending network. Section~\ref{sec:background} briefly presents
the general definition of Hawkes Processes which is then instantiated and
adapted to the contagion of gun violence in Section~\ref{sec:model}.

\subsection{Hawkes Processes}
\label{sec:background}

Hawkes processes are a class of multivariate self-exciting temporal point
processes originally introduced by Alan G. Hawkes in the early 1970s (CITE) and
have since been used to model a wide range of phenomena ranging from seismic
events to information spread in social networks to stock market trading
dynamics.

In a Hawkes process, the conditional intensity function at any given time $t$
(\emph{i.e} the instantaneous probability of occurrence of an event) can be
written as the sum of a constant exogenous intensity and endogenous
time-varying intensities for the events preceding time $t$.

Formally, for a $D$ dimensional Hawkes process, let us introduce the set of
events $\{(t_i, k_i)\}_{1\leq i \leq N}$ where $t_i$ denotes the time of event
$i$ and $k_i$ the dimension (or coordinate) on which it occurs. The conditional
intensity function is defined as follows:
\begin{equation}
    \label{eq:hawkes}
    \lambda_k(t) = \mu_k + \sum_{i=1}^N \phi_{k_i, k}(t-t_i),
    \quad 1\leq k\leq D
\end{equation}
where $M = (\mu_k)_{1\leq k\leq D}$ is the vector of exogenous intensities and
the functions $\Phi = (\phi_{i,j})_{1\leq i, j\leq D}$ is the matrix of kernel
functions (also known as exciting functions). For a pair of coordinates $(i,
j)$, $\phi_{i,j}$ models the temporal variations of the influence of coordinate
$i$ over coordinate $j$. The kernel functions are \emph{(i)} positive:
$\phi_{i,j}(t)\geq 0$ and \emph{(ii)} causal: $\phi_{i,j}(t) = 0$ whenever
$t<0$. In particular, this implies that the summation in \eqref{eq:hawkes}
is only over the indices $i$ such that $t_i< t$.

We refer the reader to (CITE) for a formal definition of the conditional
intensity function. We will simply use the following formula for the
log-likelihood of events $\mathcal{E} = \{(t_i, k_i)\}_{1\leq i\leq N}$ given
$M$ and $\Phi$ and observation period $[0, T]$:
\begin{equation}
    \label{eq:likelihood}
    \mathcal{L}(M, \Phi\given\mathcal{E}) = \sum_{i=1}^N \log\lambda_{k_i}(t_i)
    - \sum_{k=1}^D\int_{0}^T \lambda_k(t) dt
\end{equation}

\subsection{Contagion of Gun Violence as a Hawkes Process}
\label{sec:model}

We model the contagion of gun violence as a Hawkes Process by making the
following identifications: each network vertex in (\emph{i.e} each
individual) is a coordinate of the Hawkes Process and each gunshot injury is an
event of the process occurring on a coordinate of the process, the victim of
the injury.

\paragraph{Exogenous intensity.}

We assume that the exogenous intensity is the same for all the individuals in
the network. However, we attribute the regular fluctuations of the rate of
violence observed to a seasonal effect independent of peer contagion. For this
reason, we fit a time-varying function $\mu(t)$ to the data and use it for the
common exogenous intensity (see Section~\ref{sec:background}).

\paragraph{Exciting functions.} The exciting function $\phi_{u, v}(t)$ models
the effect of person $u$ on person $v$ at time $t$ and captures two common
assumptions regarding the spread of contagions.
\begin{itemize}
    \item \emph{time:} consistent with previous models used to infer the spread
        of contagions over social networks (4, 5), we assume that the impact of
        earlier infections on future events decays as the time passed since the
        original infection increases. Additionally, influence can only travel
        forward in time: an infection has no impact on those that came before
        it. As commonly done with Hawkes processes, we assume an exponential
        decay, obtaining the following formula for the temporal component of
        the exciting functions:
        \begin{displaymath}
            f_\beta(t) = \begin{cases}
                \beta e^{-\beta t} & \text{if $t\geq 0$}\\
                0 & \text{otherwise}
            \end{cases}
        \end{displaymath}
    \item \emph{network structure:} epidemiologists commonly assume that
        contagious events are localized and that the transmission probability
        increases closer to the source (CITE). In our case, we assume that
        violence is more likely to spread between people who are more closely
        linked in the network and measure the distance between individuals
        based on network topology. Based on previous studies of violence in
        social networks, we assume that infections are able to occur across
        a distance of up to three degrees of separation (6); people who are
        further away in the network have no effect on one another. For two
        vertices $u$ and $v$ whose network distance $d(u, v)$ is less than or
        equal to 3, we assume a decay of the form $\frac{\alpha}{d(u,v)^2}$.
        Hence, we obtain the following formula for the structural component:
        \begin{displaymath}
            g_\alpha(u,v) = \begin{cases}
                \frac{\alpha}{d(u,v)^2} & \text{if $d(u,v)\leq 3$}\\
                0 & \text{otherwise}
            \end{cases}
        \end{displaymath}
\end{itemize}
Finally, we combine the above two components by multiplying them to obtain the
exciting function: $\phi_{u,v}(t) = f_\beta(t)g_\alpha(u,v)$.

\subsection{Likelihood}

Using Equation~\eqref{eq:likelihood} and the model presented in
Section~\ref{sec:model}, we can now write the log-likelihood of observed
infection events $\mathcal{E} = \{(t_i, u_i)\}_{1\leq i\leq N}$ where $t_i$ is
the time of infection $i$ and $u_i$ is the vertex infected at time $t_i$. We
denote by $V$ the set of vertices in the network, and by $[0, T]$ the study
period.

Furthermore, since some individuals died during the study period, the
conditional intensity function only needs to be integrated until their time of
death in the second summand of Equation~\eqref{eq:likelihood}. Denoting by
$T_u$ the time of death of vertex $u$ ($T_u = T$ if the individual didn't die
during the study period), we obtain:
\begin{displaymath}
    \mathcal{L}(\mu, \alpha, \beta\given\mathcal{E}) = \sum_{i=1}^N \log\lambda_{u_i}(t_i)
    - \sum_{v\in V}\int_{0}^{T_v} \lambda_v(t) dt
\end{displaymath}

Using \eqref{eq:hawkes} and the explicit formula for $\phi_{u,v}(t)$, this can
be expanded to:
\begin{equation}
    \label{eq:final-likelihood}
    \begin{split}
        \mathcal{L}(\mu, \alpha, \beta\given\mathcal{E}) =&
        \sum_{i=1}^N \log\left(\mu(t_i)
    + \sum_{j:t_j< t_i} g_\alpha(u_i, u_j)\beta e^{-\beta (t_i - t_j)}\right)
    \\
    &- \sum_{v\in V}\int_{0}^{T_v} \mu(t)dt
             - \sum_{v\in V}\sum_{i: t_i < T_v} g_\alpha(u_i, v)
            \left(1-e^{-\beta(T_v-t_u)}\right)
\end{split}
\end{equation}

\begin{figure}
\centering
\includegraphics{hawkes-diagram}
\caption{Diagram of a Hawkes process. STILL NEED TO MAKE A FIGURE.}
\label{fig:hawkes-diagram}
\end{figure}

\section{Model Inference}

\subsection{Exogenous Intensity}
\label{sec:background}

Because the seasonal variations in gunshot rates (Figure~\ref{fig:background})
remain consistent throughout the study period, we assume these are inherent to
the infection process and not purely driven by noise or social contagion. We
model these seasonal variations by a periodic sinusoidal function.

Let $M(t)$ denote the number of infections occurring on day $t$ across all
individuals. We assume the following form:
\begin{displaymath}
    M(t) = A\big[1 + \rho \sin(\omega t + \phi)\big]
\end{displaymath}
Because we are modeling annual fluctuations, we know that the period is one
year, \emph{i.e.} $\omega=\frac{2\pi}{365.24}$. The remaining three parameters
($A$, $\rho$ and $\phi$) are learnt using non-linear least squares estimates
with the Gauss-Newton algorithm. This yields:
\begin{displaymath}
    M(t) = 3.78\left[1 + 0.43 \sin\left(\frac{2\pi}{365.24}
    t + 4.36\right)\right]
\end{displaymath}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{background.pdf}
\caption{Aggregated number of infections. Each blue dot represents the number
    of infections (fatal and nonfatal) that occurred on a given day. The green
curve is the sinusoidal function fit to the data.}
\label{fig:background}
\end{figure}

Because we do not know \emph{a priori} the relative importance of the exogenous
intensity and peer contagion, we only keep $\rho$, $\omega$ and $\phi$ from the
fitted parameters. In other words, we only keep the parameters characterizing
the seasonal fluctuations; the base amplitude $A$ of the exogenous intensity
will be inferred together with the exciting functions parameters.

Finally, we relate the aggregate number of infections to the node-level
exogenous intensity. By definition:
\begin{displaymath}
    M(t) = \sum_{v\in V}\int_{t-1}^t \mu(s)ds = |V|\int_{t-1}^t\mu(s)ds
\end{displaymath}
where we used that the exogenous intensity is shared across the nodes. Assuming
that that $\mu$ is approximately constant over the course of one
day\footnote{The time resolution in our dataset is the day, so we only need to
approximate $\mu$ at the day level.}, we get $M(t) = |V|\mu(t)$. Hence we
obtain the following form for the exogenous intensity:
\begin{equation}
    \mu(t) = \mu_0\left[1 + 0.43 \sin\left(\frac{2\pi}{365.24}
    t + 4.36\right)\right]
\end{equation}
where $\mu_0 = \frac{A}{|V|}$.

\subsection{Exciting Function Parameters}

Using the exogenous intensity from Section~\ref{sec:background}, the
log-likelihood now depends on three parameters $(\mu_0, \alpha, \beta)$, and
finding the maximum likelihood estimate of these parameters amounts to solving
the following optimization problem:
\begin{equation}
    \label{eq:mle}
    \argmax_{\mu_0,\alpha,\beta}\mathcal{L}(\mu_0,\alpha,\beta\given\mathcal{E})
\end{equation}
Unfortunately, the function $\mathcal{L}$ is not jointly concave in its three
arguments. We will however exploit the following fact.

\begin{proposition}
    \label{prop:concave} 
    The function $(\mu_0,\alpha)\mapsto
    \mathcal{L}(\mu_0,\alpha,\beta\given\mathcal{E})$ is concave.
\end{proposition}

\begin{proof}
    Using Equation~\ref{eq:final-likelihood}, it is easy to see that the second
    and third sums are linear in $(\mu_0, \alpha)$, hence it is sufficient to
    show that for $1\leq i\leq N$:
    \begin{displaymath}
        f(\mu_0,\alpha) = \log\left(
            \mu(t_i)
            + \sum_{j:t_j< t_i} g_\alpha(u_i, u_j)\beta e^{-\beta (t_i - t_j)}
        \right)
    \end{displaymath}
    is concave. For this, we see that the operand of the $\log$ function is
    linear in $(\mu_0, \alpha)$. By composition with $\log$ with the concave
    function $\log$, we obtain that $f$ is concave and thus conclude the proof.
\end{proof}

Numerically, we observed that $\mathcal{L}$ has many local optima, hence we
solved \eqref{eq:mle} using the following heuristic:
\begin{itemize}
    \item first, we did a brute force grid search to locate good starting
        points for the refining heuristic.
    \item second, starting from the best gridpoints obtained during the first
        step, we refined the solution by alternated minimization, that is
        alternatively:
        \begin{itemize}
            \item optimizing over $(\mu_0, \alpha)$ for a fixed value of
                $\beta$. For this, Proposition~\ref{prop:concave} allows us to
                use standard convex optimization machinery (gradient descent in
                our case) to solve this step exactly.
            \item optimizing over $\beta$ for a fixed value of $(\mu_0,
                \alpha)$. The simulated annealing was used at this step.
        \end{itemize}
\end{itemize}
Other heuristics were considered: using gradient descent as well for the
optimization over $\beta$, or using global gradient descent to optimize over
$(\mu_0,\alpha,\beta)$ at the same time. All heuristics led to the same optimal
solution, indicating that our initial grid search with fine enough to identify
good starting points. We obtained the following values of the parameters at
the optimum:
\begin{displaymath}
    \mu_0 = 1.19\cdot 10^{-5},\quad
    \alpha = 7.82\cdot 10^{-3},\quad
    \beta =  3.74\cdot 10^{-3}.
\end{displaymath}

\begin{comment}
\begin{equation}
\lambda_v(t) = 1.1845\e{-5} \left[1 + 0.43 \sin\left(\frac{2\pi}{365.24} t + 4.36\right) \right] + \sum_{u \in V} \frac{0.00317}{\text{dist}(u,v)} 0.0039 e^{-0.0039(t-t_u)}
\end{equation}
\end{comment}

\section{Inferring the Patterns of Infections}

\subsection{Methodology and Results}
\label{sec:meth}

\begin{figure}[h]
\centering
\includegraphics[width=.8\textwidth]{cascade-distribution}
\caption{The distribution of cascade sizes follows a power-law distribution.}
\label{fig:cascade-sizes}
\end{figure}

Given fitted values of the parameters of the Hawkes contagion model, it is then
possible to determine whether an infection event $(t, v)$ was primarily caused
by the exogenous process or endogenous peer contagion. For this, using formula
\eqref{eq:hawkes}, we simply compare the value of the exogenous intensity and 
the sum of endogenous exciting functions at the time $t$ of infection, and
attribute the infection event to the largest of the two quantities. In other
words, for our contagion model we compare:
\begin{enumerate}
    \item $\mu(t) = \mu(t) = \mu_0\left[1 + 0.43 \sin\left(\frac{2\pi}{365.24}
    t + 4.36\right)\right]$ and
\item  $\sum_{i:t_i< t} g_\alpha(u_i, v)\beta e^{-\beta(t-t_i)}$
\end{enumerate}

It is possible to go one step further: for an infection event attributed to
peer contagion, we can single out a single past infection event as its cause.
This is achieved by choosing the event $i$ for which the value of the summand
in 2. is maximum. That is, the event $i$ for which the value of $g_\alpha(u_i, v)\beta
e^{-\beta(t-t_i)}$ is largest.

We thus uncover the patterns of infections: each infection event is attributed
to either the exogenous intensity or a single past infection event.  Drawing an
edge between infection events $(t_1, u_1)$ and $(t_2, u_2)$ if $(t_1, u_1)$ is
the cause of $(t_2, u_2)$, we obtain the forest of infections: tree roots are
the infection events attributed to exogenous intensity, and tree nodes are
events attributed to peer contagion and are connected to their cause.  We note
that cycles are impossible since edges are directed forward in time.

Following the terminology found in the literature on information diffusion in
social networks, we call \emph{cascade} a single tree in the forest of
infections. The distribution of the cascade sizes extracted from our dataset
can be found in Figure~\ref{fig:cascade-sizes}. Consistently with previous
findings in other domains of application (CITE) this distributions follows
a power-law of exponent XX. Figure SX shows three cascades extracted from the
forest.

\begin{comment}
\subsection{Experiments with synthetic data}

\subsubsection{Generating networks}
Given that we do not know the true pattern of infection propagation in criminal networks, we first verify that our methods can accurately infer cascades in cases where we do have ground truth data. To test our methods, we generated a series of cascades on social networks. We then had our model determine optimal parameters and use these to guess how each infection spread. We show that our model is able to accurately determine the contagion’s parameters and to identify the path that infections took through social networks.

We first generated network structures using the forest fire model (CITE), which is known to capture the degree distribution and community structure observed in empirical social networks. We generated networks with 10,000 nodes and burning probabilities $p=r=0.3$.

We also simulated contagions on the co-offending network. Since we are most interested in ultimately understanding the diffusion process of violence on the co-offending network, it is important to first test our cascade inference algorithm here and ensure that we are able to accurately recover cascades. We removed all victim and demographic information, leaving just the network structure, and generated contagions. As shown below, we are able to accurately infer the process by which cascades spread on this network.

\subsubsection{Simulating contagions}

\subsubsection{Results}
\end{comment}

\subsection{Comments on Causality}

The interpretation of the cascades obtained by applying the methodology of
Section~\ref{sec:meth} requires a nuanced understanding of the notion of
causality. In particular, the meaning of assigning a single cause to an
infection needs to be qualified.

A common definition of causation is the so-called \emph{counter-factual} one
whose first formalization can be traced back to the philosopher David Hume
(CITE):
\begin{quote}
    \emph{\ldots we may define a cause to be an object followed by another
        [\ldots] where, if the first object had not been, the second never
    had existed.}
\end{quote}
The limits of this definition have been the subject of many philosophical
works, starting from David Lewis' paper \emph{Causation} (CITE). In particular,
the problem of \emph{late preemption} is made apparent by the following thought
experiment due to Ned Hall (CITE):
\begin{quote}
    \emph{Suzy and Billy, two friends, both throw rocks at a bottle.  Suzy is
        quicker, and consequently it is her rock, and not Billy’s, that breaks
        the bottle.  But Billy, though not as fast, is just as accurate: Had
        Suzy not thrown, or had her rock somehow been interrupted mid-flight,
    Billy’s rock would have broken the bottle moments later.}
\end{quote}
While Suzy's throw is obviously the cause of the bottle shattering, it fails
the counter-factual definition of causation, since had her throw not occurred
the bottle would have been shattered all the same by Billy's throw.

Without entering the philosophical debate about defining causation, we note
that, at least for the purpose of assigning responsibility, Billy and Suzy
should be considered jointly responsible for the bottle breaking: they were
both throwing rocks at it, and the fact that Suzy's rock reached the bottle
first is coincidental. Of course, if Suzy was throwing rocks twice as
frequently as Billy, without observing whose rock hit the bottle first, it
would be tempting to say that it is more likely that a rock thrown by Suzy hit
the bottle first, or that she is ``more responsible''. But it is still a fact
that at the time the bottle was broken, it was the target of the combined
throws of Billy and Suzy.

The Hawkes contagion model introduced in \eqref{eq:hawkes} can be
re-interpreted in light of this example: as they become infected, victim nodes
start throwing rocks at their associates with a frequency which decreases over
time and a node becomes infected when it is hit by a rock. Since we do not
observe whose rock hits first, the only thing we can say for certain is that at
the time it was infected, the individual was subject to the combined throws of
its neighbors who were infected in the past. This combined effect is expressed
mathematically by the sum in \eqref{eq:hawkes}. The Hawkes model also allows us
to order the nodes by the frequency at which they are throwing rocks.

It is now clear which interpretation to give to the cascades extracted in
Section~\ref{sec:meth}: it is a simplification where we assign the
responsibility of an infection to the candidate cause which was throwing rocks
with the highest frequency at the time of the infection. This simplification is
acceptable in the sense that it is more likely that the infection was the
effect of a rock coming from this cause. But even though it makes tracing the
paths of violence easier, a more conservative statement would instead consider
all the throws to be jointly responsible for the infection.
\end{document}