1 Introduction

Following the shooting at Marjory Stoneman Douglas High School in Parkland, Florida, and the subsequent activism of the students, there was a renewed interest in gun control policies in preventing mass shootings. We focus our efforts on modeling the number of mass shooting fatalities across years from 1982 to 2017.

For a deeper understanding of what, exactly, is meant by the term mass shooting, we turn to Mother Jones’ “A Guide to Mass Shootings in America” [1]. Before 2013, the FBI defined a mass shooting as a “single attack in a public place in which four or more victims were killed.” In 2013, a mandate for a federal investigation into mass shootings lowered this threshold to three or more fatalities. The data set from Mother Jones’ investigation into mass shootings [2] uses these thresholds.

Due to the importance of understanding such frequently occuring tragedies in the United States, we ask ourselves whether or not we can use population and finance models to properly model US mass shooting fatalities. Specifically, we fit three different partially observed Markov processes (POMP) models:

  1. the Ricker model
  2. a stochastic volatility model
  3. GARCH(1,1)

First, we observe some simple exploratory analysis

We quickly notice that mass shooting fatalities were relatively uncommon until after 2010, when they spiked to 71 in 2016 and 117 in 2017. We also notice that the sample autocorrelation plot indicates that the fatalities are relatively uncorrelated across time-lags. Despite the relatively uncorrelated time-lags, we proceed with POMP analysis since we desire find usefulness in models that account for potential latent variables such as fatality density and fatality volatility.

2 The Ricker Model

A common model in population biology, the Ricker model attempts to describe population growth and resource depletion. It models population growth as a function of population density, \(P_n\) at a given time. In the context of mass shootings in America, it is sensible to assume some time-varying density in the number of fatalities that drives the total number of fatalities. Thus, we formally state the model:

\[\begin{align*} Y_n | P_n &\sim \text{Poi}(\phi P_n)\\ P_{n+1} &= rP_n\exp(-P_n + \epsilon_n),\qquad \epsilon_n\sim N(0,\sigma^2). \end{align*}\]

where \(Y_n\) is the number of mass shooting fatalities in year \(n\).

First, we create likelihood slices for the \(\phi\) and \(r\) parameters, and then plot the likelihood surface. The model implementation and figures are presented below.

# Create pomp object for year mass shootings

skel = "DN = r*N*exp(-N);" # add in determnistic skeleton

# Add stochastic Ricker model 
stochStep = "
e = rnorm(0,sigma);
N = r*N*exp(-N+e);
"

# Add in the (Poisson) measurement model  
ms_rmeasure = "Fatalities = rpois(phi*N);"
ms_dmeasure = "lik = dpois(Fatalities,phi*N, FALSE);"
ms_statenames = c("N","e")
ms_paramnames = c("phi","r","sigma")
ms_initializer = "
  N = 1;
  e = 6;
"

# Build pomp object
ms_pomp = pomp(
  data = dat_year,
  times = "Year",
  t0 = 1981,
  skeleton = map(Csnippet(skel)),
  rprocess = discrete.time.sim(
    step.fun=Csnippet(stochStep),
    delta.t=1
    ),
  rmeasure = Csnippet(ms_rmeasure),
  dmeasure = Csnippet(ms_dmeasure),
  statenames = ms_statenames,
  paramnames = ms_paramnames,
  initializer = Csnippet(ms_initializer),
  cdir=getwd()
  )

## [[1]]
## NULL
## 
## [[2]]
## NULL

We observe that \(\hat{\phi} > 30\) and \(\hat{r}\approx 2\) may potentially maximize the log-likelihood for our Ricker model. We use the IF2 algorithm [3] to carry out a global maximization of our log-likelihood using randomized starting values. This allows us to obtain parameter estimates and a log-likelihood.

results_global <- data.frame(logLik=liks_global[,1],logLik_se=liks_global[,2],t(sapply(mifs_global,coef)))
summary(results_global$logLik,digits=5)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -147.9  -147.8  -147.7  -147.7  -147.6  -147.4

From the POMP diagnostics, we observe that the log-likelihood estimates, along with \(\hat{\sigma}\) and \(\hat{r}\) have converged to their maximum-likelihood estimates. However, we are less convinced that our estimate of \(\phi\) has converged. Due to our flux allocation being over-burdened, our IF2 algorithm used 2,000 particles (Np = 2000) with 200 iterations (Nmif = 200). Further analysis would involve increasing the number of particles to potentially obtain a better estimate of \(\hat{\phi}\). For our purposes, we will use \(\hat{\phi} = 40\) from viewing the likelihood surface plotted above. The convergence plots indicate \(\hat{r}\approx 1.8\) with \(\hat{\sigma}\approx 0.93\), in addition to a log-likelihood of -147.7.

Below, we simulate from the model with our obtained MLE estimates while comparing it to our observed mass shootings data. We see that the Ricker model appears to capture the sharp and volatile peaks of our data. Thus, our data on mass shootings appears to be well-modeled by the Ricker model, with the assumption of a latent time-varying density of shootings.

3 Stochastic Volatility Model

The second model we fit to our data is a stochastic volatility model from Breto (2014) [4], which models leverage at year \(n\), which here will be the correlation between mass shooting fatalities for year \(n-1\) and the increase in log volatility from year \(n-1\) to \(n\). Our primary motivation for fitting this model is the apparent volatility year-to-year in mass shooting fatalities; we hope that a successful model, such as this stohastic volatility model, can capture this. Thus, we create a POMP implementation of the following model:

\[\begin{align*} Y_n &= \exp(H_n/2)\epsilon_n,\\ H_n &= \mu_h(1-\phi) + \phi H_{n-1} + \beta_{n-1}R_n\exp(-H_{n-1}/2) + \omega_n,\\ G_n &= G_{n-1} + \nu_n, \end{align*}\]

where \(\beta_n = Y_n\sigma_{\eta}\sqrt{1-\phi^2}\), \(\{\epsilon_n\}\) is an iid \(N(0,1)\) sequence, \(\{\nu_n\}\) is an iid \(N(0,\sigma^2_\nu)\) sequence, and \(\{\omega_n\}\) is an iid \(N(0,\sigma^2_\omega)\) sequence. Here, \(H_n\) represents the log volatility of American mass shootings. This model is fit to the differenced and demeaned mass shootings data (plotted in the introduction).

The implementation of this model is below.

ms_statenames = c("H","G","Y_state")
ms_rp_names = c("sigma_nu","mu_h","phi","sigma_eta")
ms_ivp_names = c("G_0","H_0")
ms_paramnames = c(ms_rp_names,ms_ivp_names)
ms_covarnames = "covaryt"

rproc1 = "
  double beta,omega,nu;
  omega = rnorm(0,sigma_eta * sqrt( 1- phi*phi ) * sqrt(1-tanh(G)*tanh(G)));
  nu = rnorm(0, sigma_nu);
  G += nu;
  beta = Y_state * sigma_eta * sqrt( 1- phi*phi );
  H = mu_h*(1 - phi) + phi*H + beta * tanh( G ) * exp(-H/2) + omega;
"
rproc2.sim = "
  Y_state = rnorm( 0,exp(H/2) );
 "

rproc2.filt = "
  Y_state = covaryt;
 "
ms_rproc.sim = paste(rproc1,rproc2.sim)
ms_rproc.filt = paste(rproc1,rproc2.filt)

ms_initializer = "
  G = G_0;
  H = H_0;
  Y_state = rnorm( 0,exp(H/2) );
"

ms_rmeasure = "
   y=Y_state;
"

ms_dmeasure = "
   lik=dnorm(y,0,exp(H/2),give_log);
"

# Transform parameters to be defined on whole real line

ms_toEstimationScale = "
  Tsigma_eta = log(sigma_eta);
  Tsigma_nu = log(sigma_nu);
  Tphi = logit(phi);
"

ms_fromEstimationScale = "
  Tsigma_eta = exp(sigma_eta);
  Tsigma_nu = exp(sigma_nu);
  Tphi = expit(phi);
"

# Build pomp model for filtering and parameter estimation
ms.filt = pomp(data=data.frame(y=dat_demeaned,
                     time=1:length(dat_demeaned)),
              statenames=ms_statenames,
              paramnames=ms_paramnames,
              covarnames=ms_covarnames,
              times="time",
              t0=0,
              covar=data.frame(covaryt=c(0,dat_demeaned),
                     time=0:length(dat_demeaned)),
              tcovar="time",
              rmeasure=Csnippet(ms_rmeasure),
              dmeasure=Csnippet(ms_dmeasure),
              rprocess=discrete.time.sim(step.fun=Csnippet(ms_rproc.filt),delta.t=1),
              initializer=Csnippet(ms_initializer),
              toEstimationScale=Csnippet(ms_toEstimationScale), 
              fromEstimationScale=Csnippet(ms_fromEstimationScale)
)

Again, we use the IF2 algorithm to carry out a global maximization of our log-likelihood using randomized starting values. Again, we use Np = 2000 and Nmif = 200 while noting that further analysis entails using more particles and iterations. The results and diagnostics are below.

r.box = data.frame(logLik=L.box[,1],logLik_se=L.box[,2],t(sapply(if.box,coef)))
summary(r.box$logLik,digits=5)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -157.9  -157.9  -157.9  -157.9  -157.9  -157.9

We observe that the estimated log-likelihood appears to have converged, as have \(\hat{\phi}\), \(\hat{\sigma}_{\eta}\) and \(\hat{\mu}_h\). However, due to our computational limitations with our flux allocation being overburdened, the remaining parameter estimates might be unstable. Futher analysis would entail assuring that these parameters converge. Regardless, we can be confident that we have effectively maximized the log-likelihood and obtained a reliable estimator. We estimate the log-likelihood to be -157.9. With six fitted parameters, we obtain an AIC value of 327.8. We eventually compare this AIC value to a benchmark model, namely the GARCH(1,1) model.

We simulate from this stochastic volatility model using our parameter estimates and compare with the differenced/demeaned data:

We observe that the simulations clearly reflect the data we observed, lending credit to the idea that financial POMP models could be used successfully to describe problems in the United States such as mass shooting fatalities.

4 GARCH(1,1)

We fit the GARCH(1,1) model to our mass shooting data as a standard of comparison for the stochastic volatility model above. The model is outlined below:

\[\begin{align*} Y_n &= \epsilon_n\sqrt{V_n},\\ V_n &= \alpha_0 + \alpha_1Y^2_{n-1} + \beta_1V_{n-1}, \end{align*}\]

where \(Y_n\) is the differenced and demeaned data and \(\epsilon_{1:N}\) is a white noise process. After fitting the model with three parameters, we obtain the results below.

fit.garch.benchmark <- garch(dat_demeaned,grad = "numerical", trace = FALSE)
L.garch.benchmark <- logLik(fit.garch.benchmark)
L.garch.benchmark
## 'log Lik.' -153.7508 (df=3)
AIC(L.garch.benchmark)
## [1] 313.5015

The GARCH(1,1) model gives a log-likelihood of -153.8 with an AIC of 313.5. The AIC of this benchmark model is actually less than that of the more complicated stochastic volatility model, which gave an AIC of 327.8. Since we prefer models with lower AIC values, we can suggest that the GARCH(1,1) model might actually have more predictive power.

5 Conclusions

As noted in the introduction, we set out to determine if population and financial POMP models could acurately describe mass shooting fatalities in the United States. We have demonstrated that models like the Ricker model, stochastic volatility models, and the GARCH(1,1) model have applications far outside the domains in which they were created.

Simulations from the Ricker and stochastic volatility models with respective MLE estimates appear to replicate the data fairly well, while the GARCH(1,1) model may actually have more predictive power than the stochastic volatility model due to a lower AIC value.

As mentioned in the body of this report, future analysis includes ramping up computing power to obtain even better parameter estimates. In addition, a wider range of financial and population POMP models should be fit to the data to observe more fully the extent to which these models can be applied to mass shooting fatalities.

6 References

[1] Follman, Mark. “A Guide to Mass Shootings in America,” MotherJones, 20 July 2012, www.motherjones.com/politics/2012/07/mass-shootings-map/.

[2] Follman, Mark. “US Mass Shootings, 1982-2018: Data From Mother Jones’ Investigation,” MotherJones, December 2012, www.motherjones.com/politics/2012/12/mass-shootings-mother-jones-full-data/.

[3] Ionides, E. L., D. Nguyen, Y. Atchadé, S. Stoev, and A. A. King. 2015. Inference for dynamic and latent variable models via iterated, perturbed Bayes maps. Proceedings of the National Academy of Sciences of USA 112:719-724.

[4] Bretó, C. 2014. On idiosyncratic stochasticity of financial leverage effects. Statistics & Probability Letters 91:20-26.