A Dynamic Spatio-Temporal Stochastic Modeling Approach of Emergency Calls in an Urban Context

: Emergency calls are deﬁned by an ever-expanding utilisation of information and sensing technology, leading to extensive volumes of spatio-temporal high-resolution data. The spatial and temporal character of the emergency calls is leveraged by authorities to allocate resources and infrastructure for an effective response, to identify high-risk event areas, and to develop contingency strategies. In this context, the spatio-temporal analysis of emergency calls is crucial to understanding and mitigating distress situations. However, modelling and predicting crime-related emergency calls remain challenging due to their heterogeneous and dynamic nature with complex underlying processes. In this context, we propose a modelling strategy that accounts for the intrinsic complex space–time dynamics of some crime data on cities by handling complex advection, diffusion, relocation, and volatility processes. This study presents a predictive framework capable of assimilating data and providing conﬁdence estimates on the predictions. By analysing the dynamics of the weekly number of emergency calls in Valencia, Spain, for ten years (2010–2020), we aim to understand and forecast the spatio-temporal behaviour of emergency calls in an urban environment. We include putative geographical variables, as well as distances to relevant city landmarks, into the spatio-temporal point process modelling framework to measure the effect deterministic components exert on the intensity of emergency calls in Valencia. Our results show how landmarks attract or repel offenders and act as proxies to identify areas with high or low emergency calls. We are also able to estimate the weekly average growth and decay in space and time of the emergency calls. Our proposal is intended to guide mitigation strategies and policy.


Introduction
Emergency calls are considered a crucial tool to respond to incidents that require immediate attention, including accidents, wildfires, crimes, and medical emergencies.The information provided in these calls typically encompasses not only the description of the incident, but also its location and time, which are vital elements for a prompt response [1].In order to ensure an effective response, authorities analyse the spatial and temporal characteristics of emergency calls to allocate resources and infrastructure, identify areas of high risk, and formulate contingency plans.In this context, emergency calls' spatial and temporal analysis is crucial to understanding and mitigating distress situations.
The spatio-temporal analysis of emergency calls is relatively modest.Some papers employ GIS techniques to explore the spatial and temporal dynamics of emergency calls to further improve emergency services [2,3].Most works use statistical methods to forecast future events and to determine emergency call driving factors [4], spatial and temporal clusters [5][6][7][8], or approximating population sizes [9].In particular, Heaton et al. [10] and Li et al. [11] apply spatial and spatio-temporal point processes, a discipline within spatial statistics, to model the spatio-temporal characteristics of emergency calls.As emergency calls often mirror crimes, these authors adapted popular point process methodologies for crime data to analyse distress signals.For example, Li et al. [11] analysed emergency calls in Montgomery County, Pennsylvania (2016-2017) using a non-parametric spatio-temporal self-exciting point process model previously employed for crime data.The model captured the clustering features in emergency calls and distinguished the areas with intrinsically high emergency calls and those temporal intervals with higher peaks in the calls to direct emergency interventions.For a neater exposition, a summarised version of previous works can be found in Appendix A.
Spatial and spatio-temporal point processes define a suitable mathematical framework to model location-based data in various scientific disciplines, including ecology, epidemiology, and criminology.In particular, the pattern formed by the spatio-temporal coordinates of a crime in a region or a city can be represented by a stochastic point process that could be augmented with additional spatial or temporal covariate information.Point processes to analyse crime data are frequent in the literature given the natural and contextdependent tendency of crime to cluster [12].Notably, many studies use Hawkes-type point processes or log-Gaussian Cox processes (LGCPs, also called doubly stochastic Poisson processes) for modelling crime event data since both techniques account for spatio-temporal dependencies, covariate inclusion, and clustering phenomena [13,14].
Cox processes serve as suitable models for point process phenomena that are stimulated by environmental factors.However, they are not as apt for phenomena that are predominantly instigated by interactions among the points themselves.A particular property in LGCPs is that the logarithm of the intensity surface is a Gaussian process.As noted by Mohler et al. [12] and Diggle et al. [15], this produces a range of advantageous features that simplify the estimation, interpretation, and simulation of the model.Additionally, the stochastic nature of the intensity process enables the capture of spatio-temporal dependencies [16].In light of this, it can be challenging, or even unfeasible, to differentiate empirically between processes that represent stochastic, independent fluctuations in a heterogeneous environment, and those that represent stochastic interactions in a homogeneous environment [15].In the same line, Hawkes point processes, being a type of self-exciting processes, can model the space-time structure of events conditional on the history through the specification of a triggering function.
Emergency calls are defined by an ever-expanding utilisation of information and sensing technology, leading to extensive volumes of high-resolution data.These data are typically heterogeneous and dynamic, characterised by intricate underlying processes in conflict processes, such as diffusion, heterogeneous escalation, and volatility.Consequently, the temporal dynamics of crime data recorded from emergency calls can not be trivially handled by the triggering function in Hawkes processes or by classical formulations of LGCPs.We note that the latter type of processes depend on a space-time covariance structure which is difficult to handle against large datasets and complex non-separable structures.This poses both a mathematical and a computational problem.In this line, Hawkes processes cannot easily address complex dispersion processes such as advection and diffusion [17,18] and these are basic characteristics associated with emergency calls, which have to be prudently introduced into the modelling framework.To account for a system's complex temporal dynamics and to reinforce the discrete-time series definition in LGCPs, Zammit-Mangion et al. [19] introduced stochastic integro-difference equations (SIDEs) as a way to fill the existing gap in modelling complicated latent factors.
As many sorts of emergency calls (such as different types of crimes) share patterns and trends in space and time as armed conflicts, and events are generally registered in discrete-time format, we build upon the reasoning of Zammit-Mangion et al. [19] to study the spatio-temporal dynamics of several emergency calls in Valencia, Spain, for ten years (2010-2020).Our ambition is to model and predict the spatio-temporal behaviour of emergency calls in an urban environment to support criminal activity mitigation strategies and policy responses.Our strategy differs from the existing literature on LGCPs and Hawkes processes in the way we consider the transition from one time to the next one, as rather than depending on covariance structures or triggering functions, we consider integro-difference equations that are able to mimic complex latent processes, and, in turn, are able to model fractional growth or decay.This paper indeed improves on Zammit-Mangion et al. [19]'s proposal in several aspects.We consider crimes in much smaller regions, such as the street network of cities, that make the spatio-temporal interaction rather different from much larger regions.This is indeed a step forward as LGCPs are very scarce for network data.We then consider an analytical procedure for parameter selection rather than a more empirical-based approach.We finally make a finer and more friendly implementation of the methodology.Our modelling strategy is able to account for the intrinsic complex space-time dynamics of some crime data on cities by handling complex advection and diffusion processes.
The structure of this paper is as follows.Section 2 presents the data and the motivating problem and details the spatio-temporal point processes and SIDE methodology.Section 3 presents the data analysis, and Section 4 is devoted to conclusions.

Data
The dataset consists of the geocoded locations and times of the calls reported to the 112 emergency telephone number of the Valencian agency for security and emergency with headquarters in the city of Valencia (Spain).The data relate to individual crime activities from 1 January 2010 to 31 October 2020, in Valencia, Spain.Valencia, located on the southeastern coast of Spain, is one of the most important cities in the Valencian community of Spain.According to the national statistical records of Spain in 2018 (https://www.ine.es),(accessed on 30 November 2022) the city has a population of around 0.8 million.We use data from 2010 to 2019 to fit our model while holding out the first ten months of 2020 for prediction purposes and model validation.We note that we also have access to data until mid-2020.However, data from 1 March onwards cannot be used to validate the model due to, first, the quarantine, and then the severe restrictions imposed in Spain due to the COVID-19 epidemiological situation, which changed the natural effects and structure of the calls and the criminal behaviour itself.
During these ten years, 83,379 calls were registered by the emergency number in Valencia.These calls can be divided into four subgroups due to the nature and reason of the call.Among them, we have 51,533 calls due to assaults, 23,282 due to robberies to individuals in the streets, 388 due to aggression against women, and the remaining 8176 calls were due to a cause other than the previous ones but without specifying the reason for the call.
Figure 1 shows weekly data from 2010 to 2020.We note that the first weeks of 2020 show a similar trend to previous years.However, from the ninth week (first week of March), just when the mandatory quarantine for Spanish citizens was implemented, the number of calls decreased considerably to reach levels well below the rest of the years.For this reason, we work only with data until the end of February 2020 (so we only consider the first eight weeks of 2020).Figure 1 highlights a slight upward trend and similar behaviour each year.We note a rise in the intermediate weeks of the year corresponding to the summer holiday period, which is somehow expected for a touristic city.Additionally, we observe a constant peak around weeks 10-11 each year, corresponding to the weeks when the city's local festivities (well-known all over Spain) call for many tourists.The city has been divided into 81 neighbourhoods or districts, and the spatial distribution of the calls per district is given in Figure 2. We note a higher number of calls in central districts compared to others located outdoors of the city.With the idea of showing some underlying aspects of the number of calls, we now show in Figures 3-5 the weekly number of calls in six selected districts that we will use in the rest of the paper to show our analytical and prediction results.We show two neighbourhoods in the centre of the city (with a higher number of calls), Russafa and Sant Francesc, two neighbourhoods in the maritime east of the city, Cabanyal-Canyamelar and Malva-Rosa, and two neighbourhoods located on the outskirts of the city, one in the north, Benicalap, and one in the south, La Torre.Maritime districts highlight the effect of the arrival of tourists in summer with an increase in the number of calls.In addition, the central districts show a clear peak of calls in March due to the local festivities whose activities gather people in downtown Valencia.To have a graphical idea of the spatial pattern of the calls due to the selected crimes, we represent such point patterns for 2012 (see Figure 6).Here, we observe an increase in calls in March and during the summer months.In addition, we can see that the increase in March happens in central districts, while the increase in calls in the summertime is more concentrated in maritime districts.For comparative purposes and to graphically analyse the space-time interaction of the calls, we also show in Figure 7 the month of January over the ten years.This graphical output shows an increasing number of calls per year for January, giving light to such dynamics and interaction in space and time.
A final piece of information we have associated with the space-time locations of the crime-related calls is based on several covariates that are minimum distances from an event (a call) to a set of selected landmarks in the city.This is highly important as different landmarks increase or decrease the number of calls, which means criminal activity can be related to being closer or further from a particular landmark.In this line, we considered the following selected landmarks or points of interest: ATMs, banks, bars, coffee shops, industries, markets, nightclubs, police stations, pubs, restaurants, and taxi stops.Furthermore, we have measured the nearest distance from a call to such landmarks.We only show some distributions and corresponding patterns for some graphical testing for brevity.Thus, on the one hand, we see how the minimum distances to restaurants or pubs (Figure 8) are dominant in the lower values of the distribution of the distances.On the other hand, in the cases of police stations or ATMs (Figure 9), the minimum distances are not that short, indicating that crimes are happening much farther from these landmarks.Finally, in the industries or markets (Figure 10), the distances seem equal or constant, indicating these landmarks have no or little effect on crimes.

Methodology
We follow the logic of Zammit-Mangion et al. [19] to propose a suite of dynamic spatiotemporal modelling tools to identify complex underlying processes of crime data coming from calls to the emergency phone in an urban environment.Our modelling approach sits in the family of LGCP as it provides great flexibility over the more simple inhomogeneous Poisson process.Poisson point processes focus on modelling event-based data by assuming a Poisson distribution to the probability of observing a particular number of events within a defined area, D. The mean of this distribution is determined by the integral over D of an intensity function, denoted as λ(s), which is dependent on the location vector s, belonging to D. A doubly stochastic or Cox process is defined if the intensity function is assumed to be a random function.In the case of a LGCP, the logarithm of the intensity is assumed to be a Gaussian process (GP), defined as a latent structure.
Let us assume we have a discrete time series of continuous spatial LGCPs, since the temporal range is discrete.Formally, let k ∈ K, K = {1, . . ., K} denote a discrete time index set, and , a set of temporally correlated spatial Gaussian processes, each with mean µ k (s) and covariance function σ 2 k Ψ k (s, r).The intensity function of the point process is represented by a exponential function of z k (s) for each k, i.e., λ k (s) = exp(z k (s)).To mitigate prediction uncertainty, the mean function of z k (s) can be linked to explanatory variables.In this scenario, a vector of spatially referenced covariates denoted as d(s) and its corresponding regression coefficients represented by b T can be employed.As such, the intensity of the LGCP at time k is given by the exponential of the sum of the regression coefficients and the mean function, i.e., λ k (s) = exp b T d(s) + z k (s) .
In order to account for the intricate temporal dynamics of the intensity functions through z k (s), we adopt the stochastic integro-difference equation (SIDE) framework.This flexible modelling approach is capable of capturing dynamic temporal effects such as diffusion and dispersal.Specifically, the SIDE relates the spatio-temporal dependent variable z k (s) to z k+1 (s) through the following integral equation: where k M (s, r) is the mixing kernel in the integral, and e k (s) ) is an added disturbance, modelled as a Gaussian field with mean µ Q (s) and covariance function k Q (s, r), and D is the spatial domain under study.The non-linear mapping f (•) distorts the field in the sedentary stage; in the absence of a priori knowledge, the identity map f (z k (r)) = z k (r) can be adopted.

Correlation Analysis
We need to measure the correlation between the crime events within the same and across subsequent time frames.In point process statistics, these are quantified through the pair correlation function.Indeed, the pair auto-correlation function (PACF) g k,k (s, r), and the pair cross-correlation function (PCCF) g k,k+1 (s, r) quantify the probability of finding an event at location r given that an event has occurred at s within the same time frame k or at previous time frame k − 1.These functions are given by , where λ (1) ] are real and positive.The PACF determines qualitative characteristics of the events; if g k,k (s, r) = 1, no spatial pattern can be inferred from the data; g k,k (s, r) > 1 and g k,k (s, r) < 1 indicate event clustering and inhibition, respectively.
Given a spatial point pattern at time k, P k , a realisation of a particular point process model, a standard non-parametric kernel estimator of the first-order intensity λ (1) where , where w(s i , s j ) is the fraction of the circle (in two dimensions) with centre s i and radius s i − s j lying in D.
If the processes are taken to be second-order stationary also in time, to smooth out the non-parametric estimates an average over all K time steps may be taken so that gk,k , with * being the convolution operator.The above indicates that the kernel k M may be acquired by performing a deconvolution on the previous equation, and conventional image processing methods such as direct inverse filtering may be applied to achieve this.Furthermore, one can show that

Dimensionality Reduction
A computationally convenient truncated basis function representation of the spatiotemporal field is considered to develop an inferential approach.The choice of basis functions relies on the non-parametric estimation of the PACF.Specifically, we recall here two results.In the fundamental lemma of LGCPs, the log PACF equals the field auto-correlation function, The second result represents the auto-correlation theorem which indicates that the spectrum of the signal is the Fourier transform of the auto-correlation function.This connection between the frequency content of the point process and the PACF is employed to pick a suitable collection of basis functions.Having the basis functions, the kernel, the mean disturbance and the field can be decomposed as where φ(s) ∈ R n is the vector of basis functions, x k ∈ R n and ϑ ∈ R n are weights which reconstruct the spatio-temporal field and the disturbance mean, respectively.Furthermore, Σ M ∈ R n×n and Σ Q ∈ R n×n reconstruct the kernel covariance function and the disturbance covariance function, respectively.In this context, the SIDE of Equation ( 1) can be rewritten as In this equation, A(Σ I ) ∈ R n×n and w k ∈ R n represents a Gaussian coloured noise term with a mean of E[w k ] = ϑ and a covariance of cov[w k ] = Σ Q .The objective is to estimate the unknown parameters θ = ϑ, ΣM, ΣQ −1 and the states X K = x0 : K = xkk = 0 K using the data Y K = ykk = 1 K , where each y k is a set of coordinates of the logged events at the k-th time point.

Basis Function Selection
Based on the link between the Fourier transform of the PACF and the signal spectrum, we select the collection of basis functions using a frequency-based approach.The Fourier transform of the average PACF is computed, and a cut-off frequency of ν c cycles/unit is selected from it.Then, localised reconstruction kernels are placed at small regular intervals throughout the spatial domain.The centres of the basis functions {ζ i } n i=1 equate the sequence of vectors defining the regular partition of length ∆ s over the spatial domain D, so that where α 0 is an oversampling parameter.Gaussianity makes mathematical development easier while still producing flexible close forms.Thus, if basis functions are defined as Gaussian radial basis functions (GRBFs) with functional form φ(s) = exp −s 2 /2σ 2 b , their Fourier transforms are also Gaussian radial functions so that the spatial and frequency variances relate through the mappings To ensure appropriate reconstruction, the frequency and spatial range of the basis functions needs to exceed that of the field.This condition is addressed when σ ν = 1 √ 2 ν c .Given σ ν , the previously stated relation between σ ν and σ b can be used to find the width of the desired Gaussian radial basis functions (GRBFs), obtaining that σ 2 b = 2ν 2 c π 2 −1 .The resulting basis functions are a set of GRBFs with parameter σ b placed in the spatial domain D centred on the coordinates {ζ i } n i=1 .However, since the GRBFs are not of compact support, a compact GRBF (CGRBF) is used instead, which takes the form for τ > 0, and where • denotes the Euclidean distance on D. The CGRBF is similar to the usual GRBF with φ(s) = exp −τ 2 s 2 /2π but it is of compact support.The CGRBF parameter τ can be estimated having the GRBF parameter σ b , or a cutoff frequency ν c , as follows: The cutoff frequency of ν c is selected from the Fourier transform of the average PACF.The value corresponds to where the average PACF function decays abruptly towards zero.

Variational Bayesian Inference
We use the following likelihood function for inference: where each λ k (s) is approximated using the same basis representation The full posterior distribution is approximated through the variational Bayes method, and takes the form with p(•) being the variational marginals.The variational marginals inform about important properties of the crime events' progression.X K reconstructs the spatio-temporal field at every time point, ϑ shows the spatially varying escalation in events, Σ M displays the extent of the spatial dynamics, and Σ −1 Q informs about the volatility of the event occurrences.The number of unknown parameters in the reduced model scales as D(n 2 ), where n is the number of basis functions retained.
The variational Bayes marginals for the unknown states X K and parameters θ = ϑ, Σ −1 with d denoting the number of covariates, can be estimated by finding the lower on the marginal likelihood.We then have where θ /θ denotes the set of variables θ without θ and E p(•) [•] is employed to compute expectations with respect to the distribution in question.
In the upcoming section, we present the method utilised to deduce the unknown variables.It should be noted that the notation i|j denotes the estimate at time i based on the data observed until time j.For the sake of clarity, we have reformulated the model as follows: xk + 1 = xk + ϑ + wk, where wk has a zero mean.

Parameter Estimation
Starting with the state inference, the distribution p(X K ) can be computed by an ap- proximate variational Kalman smoother.Let x 0 ∼ N x0 (µ 0 , Σ 0 ).Considering the variational forward α(x k ) = p(x k | y 1:k ) and backward β(x k ) = p y k+1:K | x k messages, and using the Laplace method approximation, we can further write α(x k ) → N x k xk|k , Σ k|k and The two messages are then combined to give the smoothed estimate In relation to escalation inference, and considering the prior p(ϑ) ∼ N ϑ θp , Σ ϑ,p , its posterior p(ϑ) can be written as Considering now volatility inference, let the prior p denotes a Wishart distribution with V a positive definite, symmetric scale matrix and d degrees of freedom.The variational posterior can be then written as where the evaluation of Γ requires evaluation of the cross-covariance matrix in addition to the usual posterior covariance matrices.The computation of the cross-covariance also requires Laplace approximations.Finally, in terms of the regression parameters, under the variational Bayes approach, we let p(b) = ∏ d i=1 p(b i ), and set the prior Note, finally, that the estimation of p(X K ), p Σ −1 Q , p(b i ) requires the Laplace approx- imation.We refer to Zammit-Mangion et al. [19] in their supplementary information for further technical details.

Prediction
Assuming a linear relationship from t to t + 1, the prediction of number of events Ŷi,t+1 for i = 1, ..., L neighbourhoods and target time t + 1 immediately posterior to time t is given by Ŷi where N i,t+1 and N i,t are the estimated number of events derived from the predictive Algorithm 1, and Y i,t corresponds to the number of reported events for time t and neighbourhood i.
3: Integrate the interpolated sample over each i neighbourhood to obtain ẑk,i .4: Estimate the intensity λk,i , and average over fixed predefined intervals to obtain λt,i and λt+1,i 5: Generate two samples N i,t+1 and N i,t from Poisson random variables with intensity parameters λt,i and λt+1,i 6: Predict Ŷi,t+1 using Equation ( 4) end for

Computation of statistics
Calculate mean, median, and standard deviations out of the L estimations of Ŷi,t+1 .

Experimental Set-Up
The VB algorithm (Algorithm 2) was assumed to have converged when the change in θ and b i , i = 2, 3, . . ., 5 in subsequent iterations was less than 0.005, and when all diagonal elements in E Σ −1 Q = d V changed by less than 1%.Note that the prior scale matrix V p and the background rate b 1 arise from the observed data itself.V p was chosen such that its mean is 16I; this value is the squared reciprocal of the standard deviation of the week with the highest variance in the Levene's test for homoscedasticity.In particular, b 1 was set to −4.0, indicating the expected weekly events per year.The coefficients of the covariates b i and their variances were initialised in 0. We set the VB algorithm to run for 200 interactions; however, it usually converged between the 50th and 65th iterations.
The exploratory analysis, VB algorithm, and predictions were implemented entirely in MATLAB R2020a.We employed parallel computing, statistics and machine learning, optimisation and mapping toolboxes, and customised functions for most of the abovementioned methods.Further details about functions and parameters can be found in the code documentation at https://github.com/DavidPayares/ValenciaCallsSIDE(accessed on 17 January 2022).The VB algorithm was trained in a Windows CPU equipped with 16 GB RAM and 11th Gen Intel(R) Core(TM) i7 processor.Overall, the training time of the VB algorithm was 8 h and 32 min.Other methods' computational time was negligible.

One of the main premises of the SIDE-driven
LGCPs methodology to model the temporal dynamics is that the increments between two consecutive times are normally distributed, and thus, the system can be expressed as a Geometric Brownian motion (GBM).The GBM is redefined in terms the mean µ Q (s) and covariance function k Q (r,s) as in Equation (1).One obtains the random walk model in Equation ( 2) by decomposing the field z k+1 (s).Indeed, considering that the intensity of the LGCP at time k is given by λ k (s) = exp b T d(s) + z k (s) , we have dλ k (s) = R(s)λ k (s)dk + λ k (s)dW k (s), with the increment dW k (s) a Gaussian process with zero mean and covariance function k Q (r,s), and R(s) a spatially varying drift.
We analysed the weekly number of emergency calls in Valencia between January 2010 and December 2019.As mentioned in the data description section, the weekly behaviour is similar yearly with a general increasing trend throughout the studied period.When we examined the increments between week N and N + 1, we found that the data's 2.9% (12 weeks) were outliers.After removing these data, we confirmed that the increment rates between adjacent weeks were normally distributed.Figure 11 shows the distribution of the increments as well as their normal probability plot, suggesting normality in the temporal increments.

Basis Function Selection
In order to select a set of basis functions that describe the spatio-temporal field of the LGCP, non-parametric estimations of both the PACF and PCCF were conducted (see Section 3.1).Furthermore, the relationship between the PACF and PCCF is essential to determine the mixing kernel k M (s, r) and the noise kernel k Q (s, r) (see Algorithm 3).
Algorithm 3 Analysis for dynamic, homogeneous, isotropic spatiotemporal point processes |D| for stationary systems or simple regression where clear trends are evident.N k is the cardinality of a spatial point process P k at point k.2: Estimate ĝk,k (ν), ĝk,k+1 (v)∀k.3: Estimate k M (v) using ḡk,k (v) and ḡk,k+1 (v) 4: Estimate kQ (v) using k M and ḡk,k (v).
Figure 12a shows the non-parametric estimations of both the average of lng k,k (s, r) (PACF) and the average of lng k,k+1 (s, r) (PCCF) in terms of υ = s − r using the expressions in Section 3.1.Note that υ corresponds to approximately 270 m.Both functions are symmetric concerning zero (υ = 0).We also note that the behaviour of lng k,k (s, r) and lng k,k+1 (s, r) are almost identical; property also noticed in Zammit-Mangion et al. [19].Once lng k,k (s, r) and lng k,k+1 (s, r) have been estimated, the mixing kernel k M (s, r) and the noise kernel k Q (s, r) can be inferred using the exact inverse filter.Figure 12b shows the estimated mixing kernel of the process.This kernel should closely resemble the true underlying kernel.
Figure 12c displays a positive cross-section of lng k,k (s, r) and its corresponding confidence interval and that of the isotropic basis function selected for the modelling.φ(ν) directly comes from Equation ( 3) with the cut-off value ν c obtained from the average PACF.We chose a cut-off frequency of ν c = 0.22 cycles/units giving a basis parameter τ ≈ 1.7325 and an oversampling parameter of α 0 = 1.2 for the placement of the CGRBFs across the study area.As in Zammit-Mangion et al. [19]'s scheme, 256 basis functions were placed on a 16 × 16 grid covering Valencia.The centre of the basis functions was separated by ∆ s = 1.9 grid units.The functions were then filtered out to remove non-representative areas having sparse events.Only basis functions whose intensity was above a constant background event rate of b 1 = − 4 (approximately five reported events per year with a distance of 450 m from its centre) and whose centre does not exceed 300 m beyond Valencia's boundary were chosen.Figure 13 shows the distribution of reported emergency calls across Valencia and the basis functions selected for the study.In total, 129 basis functions meet the criteria; they cover the areas with high rates of emergency calls.Some areas in the northern and southern areas of the city do not have a basis function representation, given their low event rates.Nonetheless, we assigned to these areas the baseline background rate b 1 to ensure identifiability.

Fixed Effects
Crimes vary significantly in space and time within a region due to many demographic and economic factors, and while it is known that the combined effect of these factors favours criminal behaviour, their geographical character defines crime locations.Typically, crimes occur in business sectors within low and middle-income neighbourhoods.These sectors concentrate most of the facilities (e.g., banks, ATMs, restaurants) that compose neighbourhoods' economic activity.Offenders target victims close to locations that can maximise their profit and reduce the risk of apprehension [20,21].In this context, we have included putative geographical variables into the spatio-temporal point process modelling framework to measure the effect deterministic components exert on the intensity of emergency calls in Valencia.
We introduced ten covariates in the form of distances to relevant landmarks in Valencia.The landmarks included financial facilities such as banks and ATMs, and leisure places such as bars, pubs, restaurants, cafes, and nightclubs.We also included industrial areas, taxi stop areas, and markets.We measured the degree of correlation between the landmarks' distances and the intensity of the 112 calls.Half of the covariates displayed evidence of association with the emergency calls: distance to banks, ATMs, bars, cafes, and restaurants.
Figures 14 and 15 display the spatial distribution of distances to landmarks within Valencia and their relationship to the 112 calls' intensity values.Overall, short distances to landmarks occur primarily in the city centre and densely populated neighbourhoods; facilities, shops, and venues are located strategically to provide location advantages (e.g., access to services and amenities) for residents and tourists.Distances to restaurants are short throughout the city except in the northern and part of the southwestern neighbourhoods.
Emergency calls, which in our case reflect criminal behaviour, are known to occur frequently near facilities with features that offenders find attractive [20]; for example, places with multiple or desired targets, and our results corroborate this assumption; while the correlation between the emergency calls log intensity and the distances vary differently for each landmark, the behaviour of these relationships is relatively similar.One would expect many emergency calls in locations close to banks, ATMs, bars, cafes, and restaurants.Figure 14 shows that high intensities concentrate in a radius of approximately 1 km to 1.5 km centred in the facilities; beyond these distances, the emergency calls' intensity decreases.We notice that the standard deviations of the correlation functions (red dotted lines) widen with distances exceeding 2 km.This occurs given the sparse number of emergency calls in low crime incidence areas.In order to account for the factors contributing to the spatial intensity of emergency calls in Valencia, the distances to various amenities, including banks, ATMs, bars, industrial areas, markets, taxi stops, cafes, pubs, nightclubs, police stations, and restaurants were considered as the deterministic component of the intensity function λ k (s).This was due to the observed strong association between these distances and the average intensity of emergency calls in the region, while we introduced eleven covariates as fixed effects in the intensity function, six presented a regression coefficient of zero.Our model could identify covariates with no quantitative influence over the emergency calls intensity.Following Algorithm 2 for inference, we found the regression coefficients for only five of the distance covariates we introduced in the intensity function; the distance to banks and ATMs parameters confidence intervals were −1.9 × 10 −2 ± 1.2 × 10 −11 and −3.4 × 10 −2 ± 6.8 × 10 −11 , respectively.The results indicate that emergency calls tend to occur in proximity to economic facilities.This proximity suggests that potential perpetrators may identify victims in these areas, potentially leading to some financial gain.The regression parameters and the confidence interval for distances to bars and cafes were 0.5 × 10 −2 ± 7.3 × 10 −11 and 1.3 × 10 −2 ± 7.7 × 10 −11 , respectively.In this case, emergency calls are located far from these places.It is important to note that while these types of landmarks attract a large number of potential victims, they also increase the likelihood of apprehension by law enforcement.For example, bars and cafes often have enhanced security systems due to their elevated risk of crimes, such as robbery and assault Weisburd et al. [21].The coefficients for the distance to restaurants were found to be negative, with a magnitude of 4.0 × 10 −2 ± 6.2 × 10 −11 .This suggests that emergency calls tend to occur in proximity to restaurants.This is particularly worrying, as restaurants are often targeted by criminals for robbery, burglary, and theft, due to the accumulation of significant amounts of cash during daily operations.
These results show how landmarks, to some extent, attract or repel offenders and act as proxies to identify areas with high or low emergency calls.

Heterogeneous Growth and Decay
A spatio-temporal analysis of events is concerned with identifying high-intensity spots and their evolution over space and time, and while traditional cluster analysis excels in determining hotspots, it cannot portray the temporal characteristics that govern advection and diffusion processes.Our methodology allows us not only to locate hotspots but also to determine their behaviour over time.
Figure 16 presents the weekly average fractional growth and decay of emergency calls in Valencia.As anticipated, the majority of the city has experienced a rise in the frequency of emergency calls.However, Sau Pau, Malilla, and La Torre neighbourhoods display the highest increments over the years.For example, Sau Pau is a neighbourhood prone to robbery; by 2015, it had accumulated approximately 9% of all robberies in Valencia Las Provincias [22].Despite the fact that the central neighbourhoods of the city have the highest incidence of events, they did not exhibit a marked increase in the number of emergency calls.Conversely, areas with fewer events, such as the Quatre Carreres district, have experienced an increase in emergency call hotspots over the course of the study period.The spatial intensity of emergency calls has shown a decrease in certain areas throughout Valencia.An interesting finding is the decay in the Benicalp neighbourhood.It is considered one of the most notorious neighbourhoods in Valencia due to high poverty levels, illegal occupation, and drug traffic, and while criminal activity has grown consistently in this neighbourhood over the last few years, the number of reported emergency calls has decreased.A possible explanation is that victims or witnesses do not report crimes as they are afraid of repercussions by organised crime.Other neighbourhoods, such as La Punta and El Castellar, also display a reduction in hotspots.However, volatility (Figure 17) suggests that the data in these areas are not very reliable.

Volatility
The volatility (Σ Q ) in the SIDE-driven LGCP allows us to measure the accuracy of future intensity estimations.The lower the value in the diagonal of Σ Q , the more accurate the predictions are, and conversely.Figure 17 shows the volatility map for emergency calls in Valencia.High volatility is present in the neighbourhoods of La Punta, Benimamet, La Llum, and Castellar-L'Oliveral.In both Banimamet and La Llum, the volatility is high as the neighbourhoods reported zero emergency calls in most of the weeks of the study.This produces a volatile temporal trend fluctuating between zero and the number of logged events (e.g., La Llum reported a maximum of four events per week) that the model cannot capture, while La Punta and Castellar-L'Oliveral also report many weeks with zero observations, high volatility results from scattered events in both time and space.For example, emergency calls in La Punta are primarily located close to the boundary with adjacent neighbourhoods.The majority of the city displays low volatility values, particularly areas with generous data, such as the city centre and surrounding neighbourhoods.The volatility map suggests that we will obtain less accurate predictions as we move from the city centre towards the suburbs.

Model Fitting
We perform model parameter estimation through Bayesian inference as detailed in Section 2.2.3, and using Algorithm 2. In particular, we estimate the intensity quartiles through the posterior smoothed estimate, the smoothed covariance matrix, and the effect of the covariates.Figure 18 shows the fitted model for five different neighbourhoods in Valencia.Overall, the model fits the data well.Real values are essentially contained in the 90% confidence intervals except in weeks when the number of reported emergency calls drops down to zero, as in La Malva-Rosa neighbourhood.The model also captures the temporal trend of the events, in some cases, even when abrupt spikes occur.An example is the neighbourhood of Arrancanpins, whose events count shot up from 9 calls in week 376 to 66 in week 377.The model accurately imitates this peak.
How accurately the model fits the data varies according to the weekly changes and the overall temporal trend.The larger the shifts from week to week are, the narrower the bandwidths become.Nonetheless, generally, the model fairly reconstructs the spatial and temporal character of the reported emergency calls in Valencia.

Prediction
One of the key strengths of our methodology is its ability to make predictions.Given that we have accurately modelled the spatio-temporal dynamics of emergency calls in Valencia, we can now forecast their future behaviour.To evaluate the predictive capability of the model, we used Algorithm 1 to estimate the number of emergency calls in Valencia for the first 40 weeks of 2020 (Figure 19).We selected this period aiming to assess our model's predictive robustness.The emergency calls recorded in 2020 contrast with the training data (i.e., calls from 2010 to 2019) due to the implementation of containment measures in response to the COVID-19 pandemic.As seen in Figure 19, there was a significant decrease in emergency calls between weeks 532 and 543, which coincided with the quarantine and isolation policies imposed by the Spanish government.However, there was an increase in calls as measures were relaxed.We expect our method to be robust enough to effectively model these variations based on past data, covariates, and the interaction of space and time.To stabilise the variance, we first transformed the counts into logarithms.Table 1 displays the Pearson correlation coefficients between the predicted and actual counts and log counts.These coefficients demonstrate a strong correlation between the predicted and actual values (0.87 for counts and 0.89 for log counts), indicating the model's strong predictive ability.The scatter plot depicted in Figure 20 presents a comparison between the logarithmic median prediction of the model and the logarithmic reported cases for the year 2020.The concentration of points around the ideal prediction line demonstrates the high level of correspondence between the predicted data and the observed data.However, it should be noted that the error bars display significant uncertainty in the median estimates for some neighbourhoods, such as Castellar-L'Oliveral and Ciutat De Les Arts I de Les Ciences.This higher level of uncertainty is likely due to the presence of multiple zero observations in these areas.We cannot only predict the number of emergency calls per week but also estimate their growth and decay over space and time.Figure 21 shows the histogram of the percentage of growth and decay of emergency calls for 2020.We can see that the predictions of the growth/decay rates are remarkably accurate in La Vega Baixa, Favara, La Carrasca, Ciutat Fallera, Sani Isidre, L'Hort de Senabre, and La LLum; the estimated percentage is almost identical to the observed one.These neighbourhoods are characterised by low volatility values that range from 0.01 to 0.12.The model predicts the changes with higher uncertainty in those areas with excessive zero counts (e.g., Carpesa, Exposicio, Castellar LÓliveral, and Cami Real) and, naturally, with more considerable volatility.While the predicted counts and change rates may vary, and in some cases, differ considerably, from the ground truth values, the model is rather accurate as (i) the predicted values are always contained in the 99% confidence intervals, (ii) the predicted log counts exhibit high correlation with the actual log counts, (iii) the distributions between the observed and predicted changes are virtually analogous, and (iv) predicted values come with a measurement of uncertainty.

Comparison with Alternative Benchmark Models
In order to evaluate the predictive capability of our model, a comparison was conducted between using its logarithmic predictions and those generated by three conventional and known point process modelling methods: (i) A spatial point process with a deterministic intensity function λ k (s) = exp(b The results of this comparison provides valuable insights into the strengths and limitations of our model, as well as its overall performance in the prediction of point process data.Further details on the models are described in [15,23,24]. We assessed the four models using the mean squared prediction error (MSPE) and the Pearson correlation (ρ) between the log predicted counts and the log real counts.The results are presented in Table 2.Note that the spatio-temporal LGCP coupled with the SIDE framework displays the lowest MSPE and the highest ρ.As expected, the purely spatial point process model displayed the poorest assessment metrics, due to the absence of the temporal component in the intensity function modelling.The spatio-temporal point process model and the spatio-temporal LGCP model showed similar metrics, with the latter being a superior alternative.Despite the relatively good performance of the benchmark models, they were unable to compete with our approach.The computation of models (i) and (ii) was straightforward as did not entail the use of Monte Carlo simulations.The estimation of model (iii), on the other hand, was performed using the lgcp R package, with a processing time of 3 h and 45 min.Despite the prolonged processing time of our approach compared to other available models, the results exhibit a substantial increase in prediction accuracy.This trade-off between processing time and accuracy is justified, as the primary objective of point process models is to accurately forecasts future events in space and time.

Conclusions and Discussion
Nowadays, technology makes enormous amounts of data readily available to researchers, and being able to handle such large quantities of information provides a step forward in many societal problems.One such problem is considered here as emergency calls in an urban context.Authorities must allocate resources and infrastructure for an effective response, identify high-risk event areas, and develop contingency strategies.We highlight that such emergency calls' spatial and temporal analysis is crucial to understanding and mitigating distress situations.
We have proposed a modelling framework to handle heterogeneous, dynamic, and complex underlying processes to observe crime-related emergency calls.Our strategy can account for the intrinsic complex space-time dynamics by handling complex advection, diffusion, relocation, and volatility processes.This is shown by analysing the dynamics of some emergency calls in Valencia, Spain, for ten years (2010-2029), and we believe this is a case study that can be an excellent example for many similar emergency calls in other urban contexts.
For example, other more complex scenarios can be idealised in our framework by considering SIDE with a non-linear transformation of the random field at time k in the integrand of our SIDE equation.This would complicate the model while providing a more flexible case.A combination with some deep learning methods could be imagined here.This can be the motivation for further research.
There are, however, some limitations due to the baseline assumptions imposed in the modelling approach.Note that the SIDE-driven LGCPs methodology sits on the assumption that the increments between two consecutive times are normally distributed, and thus, the system can be expressed as a Geometric Brownian motion.This must be verified a priori, and thus we can not expect that all types of data behave this way.Another sort of computational disadvantage is that the procedure involves a large number of matrix inverse calculations that can bring trouble in cases where the determinants are close to zero.
The code used in this paper has been made publicly available on the authors' Github repository (https://github.com/DavidPayares/ValenciaCallsSIDE)(accesed date on 17 January 2022).The authors have translated the LGCP + SIDE methodology into easily understandable and executable scripts, which not only replicate the results reported in this paper but also facilitate its implementation with similar datasets.The availability of the code enhances the reproducibility and transparency of the research findings.

Figure 1 .
Figure 1.Time series by weeks of the 83,379 emergency calls in Valencia (2010-2019).

Figure 2 .
Figure 2. Map of total number of calls per district in Valencia.

Figure 3 .
Figure 3. Weekly calls distribution over time at central districts in Valencia, Russafa (top row), and Sant Francesc (bottom row).

Figure 5 .
Figure 5. Weekly calls distribution over time at Benicalap, northern Valencia (top row), and at La Torre, southern Valencia (bottom row).

Figure 6 .
Figure 6.Spatial point patterns of the emergency calls in Valencia for each month in 2012.

Figure 7 .
Figure 7. Spatial point patterns of the emergency calls for the month of January of each year between 2010 and 2019.

Figure 8 .
Figure 8. Histogram and corresponding point patterns coloured by their minimum distances to restaurants (top row) and to pubs (bottom row).

Figure 9 .
Figure 9. Histogram and corresponding point patterns coloured by their minimum distances to police stations (top row) and to ATMs (bottom row).

Figure 10 .
Figure 10.Histogram and corresponding point patterns coloured by their minimum distances to industries (top row) and to markets (bottom row).
denotes the Euclidean distance on D, c D,b (s i ) is an edge-correction factor, and k b (s) is here representing the Epanechnikov kernel.Consequently, non-parametric estimators of PACF and PCCF for v = s − r are given by ĝk,k

Figure 11 .
Figure 11.Analysis of the temporal increments.The left panel displays the temporal distribution of the emergency calls over the study period (2010-2019); the central panel shows the normallydistributed weekly increments; the right panel is the normal probability plot corroborating normality in the weekly increments.

Figure 13 .
Figure 13.Spatial locations of logged events (2010-2019), and 118 basis functions placement in the city of Valencia.

Figure 14 .
Figure 14.Maps of fixed effects and empirical relationship between distance-based fixed effects (banks) and the log spatial intensity.

Figure 15 .
Figure 15.Maps of fixed effects and empirical relationship between distance-based fixed effects (atms, bars, cafes and restaurants) and the log spatial intensity.

Figure 16 .
Figure 16.Posterior mean fractional growth (left panel) and (right panel) decay of emergency calls per week in Valencia (2010-2019).

Figure 18 .
Figure 18.Model fitting for five different neighbourhoods in Valencia.For each neighbourhood, the left panel shows the spatial distribution of the emergency calls with a buffer of 100 m, and the right panel displays the weekly counts (black) with their corresponding 90% fitted confidence intervals (green).

Figure 19 .
Figure 19.Time series of the weekly calls of the first 40 weeks of 2020 in Valencia.

Figure 20 .
Figure 20.Scatter plot and error bar plot between the log median predictions and log real values for 2020.The circles in the scatter plot represent individual neighbourhoods in Valencia, and the red line refers to the ideal prediction.The error bars represent the 99% confidence intervals for the predicted values.

Figure 21 .
Figure 21.Normalised histograms of log counts in 2020 per province as obtained from MC simulations with true change (circle in blue) and sample median (circle in red).The closer the blue circle is to the red circle, the most accurate is the median prediction.
Note that if temporally averaged PACF/PC-CFs are used, the inverse filter is given as kM

Algorithm 1
Prediction algorithm for time t + 1 Number of iterations is set to L

Monte Carlo estimation of the intensity for
iteration N = 1 to L do 1: Sample the trajectory z k through p(X K ) in t 2: Forward simulate each trajectory for t + 1 using the generative model with parameters ϑ, Σ Q and b, set to

Table 1 .
Pearson correlation coefficients between the count and log predictions of the SIDE model and the actual values for the first 40 weeks of 2020.

Table 2 .
MSPE and Pearson correlation coefficients between the log real counts and log predictions of benchmark models and our method for the first 40 weeks of 2020.