Identification of Spikes in Time Series

2.1 Summary of simulation

To compare each method’s ability to identify spikes, we devised a collection of simulation studies parameterized based on violence rates in several California cities. The cities we selected had monthly violence rates that differed in terms of mean, variance, and autocorrelation.

Each city’s simulation study included (1) simulating a series parameterized to be similar to city-level violence data, (2) adding a pre-specified magnitude and number of spikes to the series, (3) applying each detection method and calculating sensitivity (the probability of detecting a spike, given that one has occurred) and specificity (the probability of not detecting a spike, given that one has not occurred) of each method, and (4) iterating the procedure 1,000 times.

We conducted simulations with a range of spike numbers and magnitudes. The magnitudes considered were 10, 20, 30, 40 and 50 percent increases over the average rate during the study time period. The number of spikes inserted into the series were integers from 1 to 10. These variations allowed us to assess the ability of each method to detect different magnitudes of spikes and to determine whether the number of spikes present in the series influenced the performance of each method.

2.2 Data description

We selected the following cities in California for our study: Berkeley, Fresno, Oakland, Los Angeles, Richmond, Sacramento, San Diego, San Francisco, and Stockton. We selected these cities because they range in population size and their violence rates have a range of characteristics.

Interpersonal violence totals by city were created by summing the total number of deaths and injuries attributable to assault or homicide from the emergency department records and patient discharge and inpatient hospitalization records from the Office of Statewide Health Planning and Development (OSHPD) and the death records from Vital Statistics. To estimate monthly rates, we divided the number of cases in each city by the estimated number of people living in each city in each month and multiplied by 100,000. The population denominators came from the intercensal and postcensal population estimates from the U.S. Census Bureau. These data capture all assaults severe enough to require an emergency department visit or hospital stay and all homicides of California residents during 2005–2012.

2.3 Simulation details

Each city had its own simulation study. First, we fit the city’s actual monthly violence series from 2005 to 2012 (96 time points) with an ARIMA model whose parameters were selected by the Aikake Information Criteria (AIC). We then simulated from this model in order to capture general properties of the series (such as mean, variance, autocorrelation, and trend).

2.4 Fitting initial ARIMA models

An ARMA model predicts current values of the response based on past values (autoregressive (AR) parameters) and innovations or past error values (moving average (MA) parameters). If differencing is required, the ARMA model is integrated and described as an ARIMA model. The standard equation of an ARIMA model is

which is commonly expressed as

where B is the backshift operator (where Bdyt=y(t−d)), ϕ and θ are polynomials of order p and q, respectively, d is the amount of differencing, and ε is a noise process with assumed distribution N(θ, σ²).

The order of the AR portion is usually referred to as p and the order of the MA portion is usually referred to as q. An ARIMA(p,d,q) model has AR degree p, is differenced d times, and has MA degree q. We used auto.arima function from the forecast package in R, which uses AIC to select the model order (Hyndman and Khandakar 2007). The auto.arima function uses a unit root test with a null hypothesis of no unit root to pick the amount of differencing required. Once the proper level of differencing is determined, a stepwise algorithm is used to select the model order, which is described in detail in the documentation for the forecast package (Hyndman and Khandakar 2007). The model coefficients are estimated using maximum likelihood. The model parameters selected for each city are listed in Table 1.

2.5 Inserting spikes

We simulated from each city’s fitted ARIMA model and inserted spikes randomly into the series, with all 96 time points equally likely to receive a spike. The number of spikes were varied from integers 1 to 10 and the magnitude of spikes considered were 10 %, 20 %, 30 %, 40 % and 50 % of the mean violence rate Each combination of spike magnitudes and numbers were run as separate simulation studies and replicated 1,000 times.

2.6 Applying spike detection methods

A spike was considered correctly identified when a method identified a spike in violence at the time point where one was inserted by our algorithm. Any spike identified at a time point where we had not inserted a spike was considered to be incorrectly identified. We summarized the performance of each method using sensitivity and specificity. Descriptions of each method follow.

2.7 ARIMA

As described previously, an ARIMA model predicts the current value of the series using a linear combination of past values and innovations. The innovations are assumed to be N(0,σ2ˆ). When using the ARIMA model to identify violence spikes, we fit the simulated series with an ARIMA model selected via AIC using the auto.arima package (Hyndman and Khandakar 2007). We then calculated the residuals and their standard deviation. Any time point with a residual value greater than two times the standard deviation of the residuals was identified as a spike. To assess whether the choice of critical value influenced the results, we replicated the analysis identifying spikes as values greater than 2.5 times the standard deviation of the residuals.

2.8 Kalman filter and smoother

The Kalman filter is a recursive data processing algorithm most famously used in engineering problems to predict trajectories based on position and velocity (Faragher 2012; Durbin and Koopman 2012). It uses the observed values of a series to update its predictions, resulting in a weighted average of the predicted and observed values where the weights are determined by the uncertainty around each. The Kalman filter uses a state space approach to time series modeling, and attempts to model a latent state that is unobserved but for which there are recorded measurements related to the state at discrete points in time. In this way, the state space approach is similar to a hidden Markov model in which the state space of the latent variable is continuous and the latent and observed variables are assumed to be Gaussian distributed.

The state space approach requires two equations, called the state equation and the observation equation. For details of the equations and derivations in this section, see Durbin and Koopman (2012).

The equation for the state of the system is

where α_t is the state of the system at time t, T_t is the transition matrix, which applies characteristics of the system at time t – 1 to generate a prediction of the state at time t, and η_t is a vector of error terms with assumed distribution N(0,Q_t). In our example, the matrix R_t is the identity matrix.

There is also an observation equation for the system, defined as

where y_t is the measured value at time t, α_t is the true state at time t, Z_t is the matrix that maps the state to the measured value, and ε_t is a vector of measurement error terms assumed to be distributed N(0, H_t).

The following matrices must be specified or estimated:

1. T_t, the state-transition matrix, which maps the state at time t – 1 to the state at time t

2. Z_t, the observation matrix which maps the true state to the observed state at time t

3. Q_t, the covariance of the state process at time t

4. H_t, the covariance of the observation measurement at time t

The Kalman filter itself consists of the following recursion equations:

where v_t is the residual of the observed minus predicted value of the (latent) state,

where a_t is the expected value of the state given past observations,

where at|t is the expected value of the state given the past and present observations,

where P_t is the variance of the state given past observations,

where F_t is the variance of the innovations, given past observations

where a_t + 1 is the prediction for the value of the state at the next time point given past observations,

where K_t, referred to as the Kalman gain, represents the change in the estimate of the state at time t after incorporating the information from the measurement,

where Pt|t is the variance of the state given past and present observations, and

where P_t + 1 is the variance of the prediction for the state at the next time point, given past observations.

The smoothing algorithm combines the results from the Kalman filtering and does backward recursion to create new smoothed estimates for the state. This is done for each time t by combining both the estimate of the state at time t from the filtering process and the residuals at times > t. The state smoothing recursion equations are:

where r_t is a weighted combination of the innovations that occur after time t – 1 and Lt=TtKt=Tt−TtPtZtFt−1.

We used the KFAS package in R to apply Kalman filtering and smoothing to the data (Helske 2014), which uses the algorithms from Durbin and Koopman (2012). For clarity, we have tried to stay as consistent as possible with the notation from these sources.

For our example, the latent state α_t is the violence level at time t and the observed data is the measured violence rate y_t at time t. To apply Kalman filtering and smoothing, we first fit an ARIMA model to the series, which was selected using the AIC, and transformed it into a state-space model formulation. The recursive filtering and smoothing algorithms described above were applied to the data, creating smoothed predictions for the violence level at each time point. Finally, we calculated residuals and identified spikes as any value greater than two times the standard deviation of the residuals. To assess whether the choice of critical value influenced the results, we replicated the analysis identifying spikes as values greater than 2.5 times the standard deviation of the residuals.

2.9 Wavelets

Wavelets are functions that oscillate around zero and satisfy ∫−∞∞ψ(x)dx=0 (Nason 2008). The wavelet transform is similar to a Fourier transform, in that for a given function f(x), we can represent the function as a sum of orthonormal functions (Nason and Silverman 1994). However, instead of expressing the function as a sum of sine and cosine or complex exponential terms like in a Fourier series, wavelet analysis uses functions that are dilations and translations of a function called a mother wavelet ψ (Nason 2008).

For example, we can define a set of wavelets

where j and k are integers, which have the effect of dilating ψ by a factor of 2^j and translating by 2^– jk, and 2^j/2 is a normalizing constant. These wavelets form an orthonormal set, which implies that a function f(x) can be decomposed into a generalized Fourier series

where j and k are integers that index the dilation number (j) and the translation number (k) (Nason 2008). The wavelet coefficients b_j,k are the inner product of f(x) and ψ_j,k(x) such that

The wavelet transform is the process of representing the function in terms of the wavelets and their coefficients. The goal of using a wavelet transform is to calculate the coefficients that allow us to approximate the series. Thresholding the coefficients allows for different levels of smoothing. Soft thresholding shrinks the coefficients around the threshold. For a more complete description of wavelets, their properties and usage, (see Nason 2008; Nason and Silverman 1994).

We used the R package wavethresh to apply a wavelet transform to the data (Nason 2008, 2010). First a discrete wavelet transform was performed according to Mallat’s pyramidal algorithm, and the coefficients were thresholded using a soft threshold. The wavelets were reconstructed by applying the inverse discrete wavelet transform. We calculated residuals and identified spikes as values greater than two times the standard deviation of the residuals. To assess whether the choice of critical value influenced the results, we replicated the analysis identifying spikes as values greater than 2.5 times the standard deviation of the residuals.

2.10 Outlier detection

We also tested the performance of an iterative model fitting and outlier detection procedure based on (Chang, Tiao, and Chen 1988; Tsay 1988).

The outlier detection process goes as follows:

1. Fit an ARIMA model to the observed time series.

2. Calculate the residuals and the variance of the residuals.

3. Compute the likelihood ratio test statistic λ_t for an additive outlier at time t

   λt=ωtˆρtˆσˆ

   where

   ωtˆ=ρt2ˆ(1−∑i=1n−tπiˆϵt)

   and

   ρt2ˆ=(1−∑i=1n−tπi2ˆ)−1

   The π-weights are functions of the estimated coefficients of the ARIMA model, and can be expressed as π(B)=ϕ(B)θ(B). These weights can also be derived by representing the model as a recursive AR model.

4. Find the maximum in absolute value of the test statistics, and compare to a pre-specified critical value.

5. If the maximum exceeds the critical value, an outlier has been detected. Remove its effect by defining a new residual for the relevant time points and recalculate the residual variance. The adjusted residuals should have the form

   ϵ˜t=ϵˆt−ωˆπˆ(B)I(t=T)fort≥T

   where I() is the indicator function.

6. Calculate new test statistics using the modified residuals and residual variance.

7. Continue until no more outliers are identified.

8. Assume the outliers identified above are known outlier times. Estimate the model parameters and the ω at each outlier time.

9. Use the parameter estimates from the models with assumed known outlier times and begin the outlier detection process again. Continue until no outliers are found.

10. If additional outliers are detected, re-estimate the parameters, incorporating the new outliers into the known outlier times.

11. Once no more outliers are identified, the procedure is complete. The locations of outliers (if any) in the time series have been identified and model parameters that exclude the effects of outliers have been estimated.

We used the detectAO function in the TSA R package to implement this outlier detection procedure (Chan and Ripley 2012).

2.11 Hypothesized performance across methods

The methods we selected have different strengths and are thus likely to perform best under different scenarios. We expect ARIMA modeling to perform best when there are few spikes and the data have simple parameterizations. When there are many spikes, the coefficients fit by the ARIMA model may be biased due to the spikes, and spike detection may suffer. While the Kalman filter and smoother assumes linear equations and Gaussian errors like the ARIMA models, the Kalman filter is more adaptive than ARIMA because of the extra uncertainty it incorporates in the measurement equation and the updating step. Therefore, we expect the Kalman filter to out-perform the ARIMA model in situations of high variance. However, because of the updating step, we expect the Kalman filter to have worse specificity as it may incorrectly adjust toward lower spike levels, obscuring their effects. In contrast to the ARIMA and Kalman methods, wavelets are useful when data have localized patterns, non-linearities, and discontinuities. Therefore, we expect these methods to perform well when the series are not well characterized by ARIMA models or when there are many spikes in a series. Outlier detection methods may perform similarly to the ARIMA and Kalman methods in situations with large shocks, but may have better performance with small shocks. The iterative procedure may capture small spikes better than other methods because they remove large spikes before searching for smaller magnitude spikes.