Simultaneous multivariate Hawkes-type point processes and their application to financial markets

In economic and financial time series we sometimes observe sudden and large price jumps. Although these events are relatively rare, they have significant impacts on not only a given financial market but also several different markets and wider macro economies. Using simultaneous Hawkes-type multivariate point process (SHPP) models, it is possible to analyze the causal effects of large events in the sense of the Granger-non-causality (GNC) from one market to other markets as well as the instantaneous Granger-non-causality (IGNC). We investigate the financial market of Tokyo and other major markets, and apply GNC tests to investigate the interdependence of large events among markets. Several important empirical findings emerge among financial markets and wider macro economies.


Introduction
In economic and financial time series, we sometimes observe sudden and large price jumps. Although they are relatively rare events, when they occur they often have significant impact on not only a single financial market but also several different markets and wider macro-economies. Several recent notable events in European and Asian countries with large jumps include the global 2008 crisis.
The standard econometric method for investigating economic and financial time series has been the statistical analysis of discrete time series. In statistical time series analysis, we often assume that the observed time series data are equally spaced realizations of stochastic processes and that the state space is d (d ≥ 1) in multivariate cases. Many statistical procedures in discrete time series analysis have been developed and applied to economic and financial time series in recent decades. When we do not observe events frequently, however, traditional discrete time series modeling with continuous state space may present important limitations. For instance, it may be difficult to distinguish major contagious events from small contagious events among different financial markets across international borders.
In this paper, we propose an alternative way of investigating economic and financial events with time series data in macro-economies, i.e., the statistical analysis of marked point processes to identify and explore multivariate time series events. Although this is not a standard approach in time series econometrics, there have been applications of this methodology in statistical seismology [see Ogata (1978Ogata ( , 2015 and the related literature, for instance]. We will show that this approach is a useful alternative way of investigating multivariate economic and financial markets to shed new light on issues that have hitherto sometimes neglected. In particular, we propose using simultaneous Hawkes-type multivariate point process models and their applications in this study. We argue that using the simultaneous multivariate Hawkes-type point process (SHPP) models, which are a new class of multivariate point processes, it is possible to investigate the causal effects of sudden and large jumps with their magnitude. We develop a new way to measure the Granger non-causality (GNC) and instantaneous Granger non-causality (IGNC) through stochastic intensity modeling of point processes.
In econometric time series analysis, the concept of Granger causality (after Granger 1969) has become an important, well-established tool for investigating relationships among multivariate time series variables (See Hosoya et al. 2017 for the recent developments). In the econometric literature, Florens and Fougere (1996) investigated several Granger causality concepts in the framework of continuous time stochastic processes, but their formulation of the problem was incomplete because they excluded the possibility of co-jumps therein, which means the simultaneous jumps that can be observed in multivariate times series data are excluded from the outset. The problem of co-jumps is important because we often analyze economic time series data in discrete time (monthly, weekly, daily, hourly, and/or minute), but continuous stochastic processes are also salient recently in financial econometrics. We need to coherently unify discrete time series analyses and continuous stochastic processes. In this paper, we investigate the possible use of co-jumps in a systematic way and develop new GNC and IGNC tests, which may provide important insight for advancing the development and application of econometric time series modeling.
Previously, Kunitomo et al. (2017) have used the traditional multivariate Hawkestype point process (THPP) models without co-jumps and the SHPP models are the extension of their models (that is, the THPP models are special cases of SHPP models). There are important cases when we need the SHPP models as we will illustrate in Sect. 6. In statistical seismology, researcher often uses the earthquake data of large magnitudes (greater than 3, say) and the time scale of measurement is short due to physical laws. In financial markets, however, the impacts of shocks occur in actual trades of financial commodities and digestion of bad or good news among market participants often needs time (1, 2 h, a day and days). Therefore, it may not be fruitful to use the method developed for seismology to financial problem mechanically and we need to consider the issues of time scales, discretization of statistical models and their measurements carefully.
Several recent studies in financial econometrics have utilized point processes and conditional intensity modeling (Ait-Sahalia and Jacod 2014; Ait-Sahalia et al. 2015;Embrechts et al. 2011;Grothe et al. 2014;and others). In a survey of these and other works, Bacry et al. (2015) noted that the focus therein is mostly on studies of financial micro-market structures. The approach developed in this paper is related to these works, but the main purpose is quite different because we develop a new point process approach to assess the relationships among different (international financial) markets. In this respect, there have also been studies on international linkage in financial markets (e.g., Hamao et al. 1990), but our proposed methodology is notably different because those studies utilized standard discrete time series modeling. In terms of empirical examples, we will investigate the interactions among Tokyo-New York, Tokyo-London, and Tokyo-Hong Kong financial markets and apply the GNC tests developed herein to those contexts. This yields several important findings among major financial markets.
The remainder of the paper is organized as follows. In Sect. 2 we present a general formulation of simultaneous multivariate Hawkes-type point process (SHPP) models. In Sect. 3, we describe the estimation method and develop non-causality tests in the sense of Granger (1969). In Sect. 4, we explore simulation results and the empirical applications are offered in Sect. 5. Concluding remarks are presented in Sect. 6 and the mathematical details are provided in the Appendix.

Simultaneous Hawkes-type point processes
We divide the observation period [0, T] into discrete periods I n i = (t n i−1 , t n i ] (i = 1, … , n) and set the initial time as t n 0 = 0 . We may interpret I n i as the i-th day, but it is possible to use higher resolution of observation periods (e.g. hourly or per minute) and we allow the irregularly-spaced time series modeling in principle.
Let the observable price processes of Itô-semimartingale with the state space of d (see Ikeda and Watanabe 1989) be P j (t) (j = 1, … , d; t n i−1 < t ≤ t n i , i = 1, … , n) , and in s ∈ I n i we denote the (negative) log-return of prices X n j (s) (s ∈ I n i ) as Let the first stopping time when X n j (s) exceeds the threshold u j (> 0) in s ∈ I n i be n (i, j, 1) . When n (i, j, 1) < t n i , we re-define the return process X n j (s) = − log[P j (s)∕P j ( n (i, j, 1))] (s ∈ I n i , s ≥ n (i, j, 1)) and let the first stopping time when X n j (s) exceeds the threshold u j crossing from below in s ∈ ( n (i, j, 1), t n i ] be n (i, j, 2) . In this way, we define the sequence of n (i, j, k) (k ≥ 1) successively.
Let also a sequence of numbers of jumps in an interval be J j (i) = #{k ∶ n (i, j, k) ∈ I n i , k ≥ 1} and then define where N j (I n l , t, u j ) is the number of counts that the resulting return process X n j (s) (s ≤ t ) exceeds u j as the threshold crossing from below in s ∈ I n l and s ≤ t ∈ I n i .
The stochastic process N n * j varies at most one in every discrete observation period.
This formulation of normalized counting in each intervals I n i (i = 1, … , n) allows us to measure the market impacts of financial price jumps in discrete intervals while we can use the standard statistical intensity modeling. We note that J j (i) are countable in [0, T] (T is finite, n is sufficiently large and u j > 0 are fixed constants) and the number of instantaneous jumps (i.e. |P j (s) − P j (s−)| > u j > 0 ) is finite for the price processes of Itô-semimartingales. In financial risk management and regulation, the sudden and/or large downward movements of financial commodity prices in short time intervals are the most important subject because of their negative consequences to financial markets and economies. 1 In the following analysis, we consider the situation that there is a common threshold value u for u j (j = 1, … , d) and the observations of the counting processes are available at day-start, day-minimum and day-end in each interval. In principle, however, we can allow the irregularly-spaced time series modeling, but we need an explicit discretization of observations. For these intervals of observation, we consider the point processes, N n * j (t, u) (j = 1, … , d), which are simple. They satisfy the standard condition for point processes that as Δt → 0 , we have the conditions where  n t− is the −field generated by the latest information before t. The (conditional) intensity functions are given by where we use the notation  n t− as the latest information before t because we discretize the counting process and we have discrete observations. For convenience, we denote  n t− as  t in the following analysis whenever such notation allows. For expository purposes, in the following analysis, we interpret the increments of N n * j (s, u k ) as if jumps of the counting process occur at t n i , the end of each interval I n i and we have set the threshold u j = u (j = 1, … , d) . When we consider the situation when the interval length goes to zero, i.e., Δ n t = max i=1,…,n |t n i − t n i−1 | ⟶ 0 as n ⟶ ∞ for a fixed T, the counting process, which is a simple point process, N n * j (s, u) weakly converges to N * j (s, u) . The resulting counting process can be interpreted as a limiting continuous-time stochastic process in high-frequency asymptotics, which is not a diffusion type but a pure jump process (see Ikeda and Watanabe 1989; Ait-Sahalia and Jacod 2014; Kurisu 2018, for example).
Next, for u j = u (j = 1, … , d) we define the point processes N n * jk (s, u) by the number of stopping times that X n j (s) exceeds u (j = 1, … , d) for a particular j, X n k (s) exceed u (k = 1, … , d;k ≠ j) for another k, and other X n l (s) (l ≠ j, k) do not exceed u crossing from below by time s in the interval I n i . By this construction, we can introduce the point processes N n * jk (t, u) with co-jumps of N j and N k by where n * jk (t, u) are the conditional intensity functions of co-jumps.
Then, when we have co-jumps of two point processes, we can define the point processes and the corresponding conditional intensity functions are given by The resulting point processes can be interpreted as the marginal point process for the j-th component of the vector point process n (s, u) with d dimension.
By extending this formulation to more complex co-jumps, in general we define where the index set J j = {j 1 , … , j l } ∈ {1, … , d} is a subset of (1, … , d) . The index sets are defined as s, u). We use the self-exciting form of conditional intensity functions J n * j (⋅) for co-jumps as n * jk (t, x| n t− ) in the same way and the marginal conditional intensity function for the j− th components as There is a one-to-one transformation between N n j (s, u) for j = 1, … , p and N n * j 1 ,…,j k (s, u) for 1 ≤ k ≤ d and between n j (t, u) and n * j 1 ,…,j k (t, u) for p = 2 d − 1.
The self-exciting Hawkes-type conditional intensity functions for the marked point processes are given by for j = 1, … , p , where N * n J i (ds × dx) are the marked point processes, j,0 are the initial intensities, g i (t − s) = e − i (t−s) are the damping functions, and C(X) = (c ji (x)) are the impact functions.
Since we are interested in sudden and large jumps of the underlying price processes (in the sense of negative returns), it is important to use the probability functions of the return process in the tail areas. In this respect, many empirical studies on the stock markets found that stock returns exhibit the non-Gaussianity and thick tails. Hence it may be appropriate to use the generalized Pareto distributions (GPD) as tail probability functions for x > u > 0 (j = 1, … , d) as (See Resnick 2007 for the details of GPD in statistical extreme value theory (SEVT)).
Herein, we assume that given the return at s X n j (s) (j = 1, … , d) , the conditional density functions are given by In terms of impact functions and the intensity function of co-jumps, there can be many possible specifications. In our empirical study, we mainly investigate the form and where a ij and c are some constants.
In particular when p = d and c ij = (i, j) (indicator functions), they correspond to the traditional multivariate Hawkes-type (THPP) processes, which are simple point processes without co-jumps.
Let p × 1 vector point process n (t, u) be partitioned as (d + (p − d)) × 1 processes as vector of co-jump point processes. The corresponding conditional intensity functions as and p × p matrices We use notation such that n 1 (t, u) is the vector process of conditional intensities of marginal jumps, diag(⋅) for diagonal matrices and we often omit n for n J i (s) (i = 1, … , p) and N n J i whenever their meanings are clear in the following analysis.
Next, we rewrite (6) and (7) as and where 1 is a d × p matrix as We call the above Hawkes-type conditional intensity models as the simultaneous multivariate Hawkes-type point process (SHPP) models because the resulting marked point processes are not necessarily simple. 2 The classical Hawkes-type point processes have been useful in applications because they are simple point processes. However, they exclude the possibility of simultaneous jumps or co-jumps and they are not fit for our purposes here. The foregoing constructions of marked point processes can be regarded as an extension of Solo (2007).

Stationarity of Hawkes-type processes
In our applications, we use stationary self-exciting Hawkes-type (marked) point processes. We take the expectation of the intensity function of (11) and (12)

The stationarity implies = [ ( (s−)] and we can use the identity relation
Then by a direct calculation, we have a set of differential equations provided the initial condition (0) and the non-degeneracy condition | | ≠ 0.
We need a condition for the convergence of (t) as t → ∞ . Therefore, the condition for the existence of stationary point processes is that the spectral radius

Applying the Bartlett spectrum
Hawkes (1971) introduced the spectral density for the stationary vector point process (t) = (N i (t)) , which was originally developed by Bartlett (1963) without cojumps ( p = d ); it is defined for the conditional intensity vector in the form of Then, when p = d (there are no co-jumps), the Bartlett spectral matrix for frequency (∈ ) is given by To permit co-jumps, Bartlett spectral matrix for the d-dimensional marginal point process vector can be defined by Then we define the relative power contribution (RPC) of the marginal spectral density function g ii ( ) (i = 1, … , d) where the frequency can be defined using the joint spectral density matrix ( ) . The (i,i)-component of ( ) can be represented as and where a ij ( ) (i = 1, … , d; j = 1, … , p) are the functions of complex variables. In addition, the instantaneous RPC ( IRPC j→i ) can be defined by In this way, we can measure the RPCs for any frequency , which corresponds to the Granger-causality measures in the frequency domain. One important aspect of the above formulation is that we have a natural definition of instantaneous Granger causality in the frequency domain, that is different from discrete time series modeling.

Conditional probability prediction
An important application of conditional intensity modeling involves assessing the conditional probability of rare events in the future from past observations. Let (j) (j = 1, … , d) be the first arrival time of an event in the j− th variable. Then we can write the probability of the random variable (j) as where  N T is the −field of information available at time T < T ′ and n j (t, u| N T ) is the conditional intensity of the j-th variable. Kunitomo et al. (2017) conducted some experiments suggesting that useful information on the conditional probability of future events can be extracted from past observations. For instance, they provided an important example vis-à -vis the conditional probability prediction of the Lehman Shock occurred in 2008 as a global crisis given past information available before that event. This illustrates the potential value of our approach.

Likelihood function
When the point process is simple and there is no co-jump, the log-likelihood function of the (d-dimensional) multivariate point process is known (see Daley 2003;Kunitomo et al. 2017) and it is given by The log-likelihood function of the marked multivariate point process with the den- 1 3 and the density function for the tail probability is given by Then we can apply the maximum likelihood (ML) method to L 1T and L 2T separately. In this formulation we use the GPD for the marginal distribution of return process.
When co-jumps are permitted, the log-likelihood function of the (d-dimensional) marginal point process is not as per the above form; instead, it should be given by where and L * 2T = L 2T . In our applications, we are principally concerned with the case in which d = 2 ; thus, there is only one extra term in the likelihood function because p = 2 d − 1.
We assume the stationarity condition (17) and the existence of second order moments of ( ) = c ij ( (s)) in the statistical inference of Hawkes-type point processes without and with co-jumps. Further, we take ( ) as the stationary conditional intensity and some q × p predictable processes (t) having second order moments. (Here q ≥ 1 and we utilize the notation T (t) in Appendix, for instance.) Then, because of the resulting martingale property given the information available at each time, it is straightforward to confirm the asymptotic properties as we have and as T → ∞.
For the one-dimensional point processes with the stationary intensity function with p = q = 1 , Ogata (1978) gave a set of sufficient conditions for the consistency and asymptotic normality of the ML estimation. His derivations are based on a martingale central limit theorem (MCLT), and it is straightforward to extend his arguments to the multi-dimensional case. For the sake of completeness, we provide details of our approach based on a new MCLT in the Appendix, which may be more general than the standard literature. In the next subsection, we develop new non-causality tests in the sense of Granger, which are explored in the context of our empirical applications.

Non-causality tests
We develop and use novel GNC tests based on the likelihood ratio principle for the Hawkes-type point processes. In particular, our results in this subsection, whose derivations are given in the Appendix, include not only the multivariate extension of existing results, but also cases in which the resultant limiting Fisher information matrix can be random variables. We first state our results for the case of no co-jumps under a set of regularity conditions, which will be extended to the more general case. The proof is lengthy, but often along the standard line of asymptotic arguments and we only give its outline in Appendix.
Theorem 1 Let the log-likelihood function of the Hawkes-type point processes with true parameters be L T ( 0 ) in (26) and (27), the log-likelihood function with the ML estimator ̂M L be L T (̂M L ) under ∈ Θ and the log-likelihood function with the restricted maximum likelihood estimator ̂R ML be L T (̂R ML ) under ∈ Θ 1 ( Θ 1 ⊂ Θ ). We assume the sufficient condition for stationarity, the existence of the second-order moment condition of ( ) , and we assume that the parameter spaces ∈ Θ in r and ∈ Θ 1 in r 1 (0 ≤ r 1 < r) are compact sets. Under a set of regularity conditions (see Theorem A-3 in the Appendix), as T → ∞, where r − r 1 is the number of restrictions of = ( k ) and 2 (r − r 1 ) is the 2 −random variable with r − r 1 degrees of freedom.
The details of a set of regularity conditions are discussed in the Appendix. When co-jumps are permitted in the Hawkes-type processes, we cannot apply Theorem 1, but it is important to obtain the corresponding results in such cases for econometric applications. When we use discrete versions of point processes, which would be often the case in econometric applications, we need to consider the existence of co-jumps. We then develop non-causality tests based on the likelihood ratio principle. In this respect, note that in our setting discussed in Sect. 2, although we permit co-jumps, it is possible to apply the martingale central limit (MCLT) theorem for point processes. Our results, in consideration of co-jumps, are an extension of Theorem 1. The proof is lengthy, but often along the standard line of asymptotic arguments and we only give its outline in Appendix.
Theorem 2 Let the log-likelihood function of the Hawkes-type point processes with true parameters be L * T ( 0 ) in (29), the log-likelihood function with the ML estimator ̂M L be L * T (̂M L ) under ∈ Θ and the log-likelihood function with the restricted ML estimator ̂R ML be L * T (̂R ML ) under ∈ Θ 1 ( Θ 1 ⊂ Θ ). We assume sufficient conditions for stationarity, and the existence of the second-order moment condition of ( ), and we assume that the parameter spaces Θ ∈ in r and Θ 1 ∈ in r 1 (0 ≤ r 1 < r)

are compact sets. Under a set of regularity conditions (see Theorem A-3 in Appendix), as T → ∞,
where r − r 1 is the number of restrictions of = ( k ) and 2 (r − r 1 ) is the 2 −random variable with r − r 1 degrees of freedom.

Simulations
To examine the relevance of the estimation and testing procedure proposed in this paper, a set of simulations are executed. The model used in these simulations is a simultaneous Hawkes-type model with two dimension and the intensity functions are given by where n 1,0 > 0, n 2,0 > 0, n 12,0 > 0 and > 0. We first generate stock price returns using the GPD as marginal and the twodimensional Gaussian copura. Then we employ ML method to obtain estimates of the underlying parameters. We provide a set of visualization (Figs. 1, 2, 3 and 4) to illustrate the key results on the finite sample distributions of the ML estimator. All histograms are standardized for the comparison of the standard normal distributions as where = ( i ) is a vector of parameters and ̂ is the ML estimator.
In our numerical evaluations, the values of estimate sometimes hit the boundaries of the non-negativity of intensity functions with finite samples, resulting in instabilities. To mitigate against this, we thus set non-negativity restrictions on parameters in our simulations. The ensuring results are reasonable, but sometimes we observe that the ML estimators of coefficients exhibit biases, although they are not very large (Fig. 2 is a typical example of this). The sample size in (34) 1∕2 n (̂− ),  our experiments was around 1000 because it may be similar to the data size in the empirical examples and it seems that we need large number of data for reducing these biases. We summarize the configuration of our numerical experiments: the number of replication of simulations as 100, and for GPD( j , j ) we set ( 1 , 1 ) = (0.007, 0.22) , and ( 2 , 2 ) = (0.008, 0.15) . These numerical values give reasonable results, and they are based on the preliminary estimates from our empirical studies. Among many simulations we illustrate key results in Table 1 and Figures. Note that because we have taken * 12 = 0 , we have a sampling distribution around zero, and the resulting estimate is not significant as illustrated in Fig. 1. Other estimates of ij , which are around their true values, take reasonable values on average in the sense that they are not significantly different from the true values, and the sampling distributions are illustrated in Figs. 2, 3 and 4. We observe some positive biases on the estimates of ij and negative biases on the estimates of initial intensities, which may be due to the results of the non-negative constraints of the parameter restrictions and the number of sample size we employed. We have imposed the non-negativities of the intensities of variables directly in the ML computation.
In the ML estimation, there can be some effects of initial conditions and we have investigated this problem in the SHPP models, where such sensitivity is also apparent, albeit minor in overall simulations.
We also use the 2 -distributions as the limiting distributions of the likelihood ratio statistics for hypothesis testing in our empirical study. We confirm that the 2 -approximations with finite samples are basically appropriate.

Empirical applications
In this section, we report the empirical results on two empirical examples using the SHPP and THPP models. The first concerns the three major stock markets, namely, Tokyo, New York, and London. Since time differences exist when each market is open and closed, it is reasonable to assume that there are no co-jumps. In terms of the second example, we focus on analyzing the simultaneous interaction among Tokyo and Hong Kong financial markets. In this latter case, since the time zone differences are small (just a 1 h difference) compared to the first empirical example, it may be natural to use SHPP, which is the extended Hawkes-type point process model with co-jumps. Because of the limitations of data available to us, we have ignored the possibilities of crossing the threshold from below except the first one in a day, or between the day-start to day-minimum.
In the first example, daily data of day-start to day-minimum data are employed covering Nikkei225, S&P500 and FTSE100 during January 2, 1990-August 26, 2015. We choose u = 2% based on the earlier study of Kunitomo et al. (2017), which used the formulation of discrete process of returns and analyzed daily data of day-start to day-end for this case. Their empirical results were quite similar to those in the following analysis, but the numerical values are different. All computations were carried out by the original programs written in R. Example 2, which concerns the Tokyo and Hong-Kong markets, is entirely new and is the principal driver of the SHPP models developed in this study. We will report the results for Example 2 using this type of data. Nonetheless, we have done robustness check of our results on the estimation of conditional intensity modeling and non-causality tests. We omit reporting the details of some results using the day-start to day-end data because they are basically quite similar.

Example 1 (Tokyo-NY-London)
We first maximize the likelihood L 2T to estimate the marginal distributions of financial market returns. As shown in Table 2, we confirmed that the marginal distributions of market returns (i.e., log-returns) have thicker tails than the normal distribution. It is because the estimates of i (i = 1, 2, 3) are positive and it is appropriate to use the GDP in our estimation procedure. The result means that the Frechet-type tail distribution is appropriate since the domain of attraction of Gaussian distribution is Gumbel (see Embrechts et al. 1998 Chapter 3 for instance). The standard deviations (SD) (or the standard errors of the estimates) in Tables are estimated by the numerical evaluation of the Fisher information matrix.
For the estimated models with two dimensions ( d = p = 2 ), we take the impact functions c(x) as Model 1 ( c(x) = 1, Model 2 ( c(x) = x, and Model 3 c(x) = x c (0 < c < 1) . The estimated values of the log-likelihood and Akaike information criterion (AIC) are those with the marginal distribution L 1T . The full likelihood can be calculated using L 1T and L 2T . The standard deviations (SD) or standard errors of the estimated coefficients are also evaluated numerically using the inverse of the estimated Fisher information matrix.

Model 1
We estimated the intensity function as  Since the ML estimates can be numerically sometimes unstable without any restrictions on the parameter space, we have set restrictions that the discounted parameters ij (i, j = 1, 2) have the same value in the following estimation. The estimation results for Case 1 are presented in Tables 3, 4. In Table 3, N 1 of Model 1 corresponds to Tokyo and N 2 corresponds to New York in Tokyo-New York markets. In terms of Tokyo-London, in Table 4 N 1 of Model 1 corresponds to Tokyo while N 2 corresponds to London.
The most important finding here (and in Tables 5, 6 below), is that the coefficient 12 is statistically significant while the coefficient 21 is not statistically significant. This represents a kind of non-causality test, but we will discuss this more formally below. We found reasonable values for other parameters in their magnitudes and significance, and they are significant for both Tokyo-New York and Tokyo-London markets.

Model 2
We estimated the intensity function as and the estimation results are presented in Tables 5, 6. In the present case, we have similar values for the estimated coefficients as Case 1 except 21 . The significance of coefficient is more pronounced here compared to Case 1, which corresponds to the likelihood values and their AIC.

Model 3
We estimated the intensity function as 0  Although we used the ML estimation in this case, the estimates of ML are often unstable numerically. In particular, we often found numerical difficulty to calculate the standard errors of estimates (Tables 7, 8). It was probably because the optimization computations with R often were unstable without any restrictions on the parameter space due to the near-singularity of the estimated Fisher information. Then we have tried to set some restrictions that the discounted parameters j (j = 1, 2) have the same value and we set c 11 = c 12 , c 21 = c 22 for instance. The results of estimation with this restriction have been given in Kunitomo et al. (2017) with the datasets of day-start to day-end. Here we report the estimation results of Model 3 with further restriction as c = c 11 = c 12 = c 21 = c 22 . Overall, the results suggest that Models 2 and 3 are better than Model 1. In addition, according to AIC, Model 2 is better than Model 3 mainly because the latter is over-parametrized for Tokyo-New York markets. Hence, we adopted Model 2 in the following non-causality tests.

Non-causality tests
In applying the GNC test procedure, we set the impact function as c(x) = x . We report our empirical results for the hypothesis H 0 ∶ ij = 0 using the likelihood ratio test (LRT) statistics based on the Tokyo-New York data. For the null-hypothesis H 0 ∶ 21 = 0, LRT statistic is 2 × (−3562.015 + 3562.017) ∼ 0 , and we could not reject the null-hypothesis. (The upper 95% critical point of 2 (1) is 3.481 in Table 5.) This means that changes of the Japanese financial market have little impact on the U.S. financial market.
For testing the null-hypothesis H 0 ∶ 12 = 0, LRT statistic based on the Tokyo-New York data was 2 × (−3562.017 + 3572.843) = 21.652, and the nullhypothesis was rejected. Thus, there is a significant effect from U.S. financial markets to Tokyo financial market (see Table 5).
Similarly, in Tokyo-London markets, for the null-hypothesis H 0 ∶ 21 = 0, LRT statistic was 2 × (−3660.215 + 3660.215) ∼ 0.0; that the null-hypothesis was not rejected. Thus, knock-on effects from Tokyo to London financial market are rather limited.
For the null-hypothesis H 0 ∶ 12 = 0, LRT statistic based on the Tokyo-London data was 2 × (−3660.215 + 3665.593) ∼ 10.756 , and the null-hypothesis was rejected. This means that the London market affects Tokyo market (see Tables 4 and  6).
To summarize our findings among three major financial markets, the effects of the Japanese market on the U.S. and London are rather limited, while we found significant effects of both of these markets on the Tokyo market. This finding agrees with several empirical findings obtained using different statistical methods as explained by Kunitomo et al. (2017).

Example 2: Tokyo-Hong Kong markets
For the second example, we have used daily data of day-start to day-minimum from the Nikkei-225 and the Hansen Index of Hong-Kong during January 2, 1990-August 26, 2015, which is the same sample period as Example 1. Since the trading periods in these two financial markets are quite similar, we expected simultaneous movements in the two markets. For the estimated models with two dimensions ( d = 2, p = 3 ), we take the impact functions c(x) as Model 1 c(x) = 1 and Model 2 c(x) = x . Because there can be many additional parameters in Model 3, which has the general form of impact functions, the estimated results are often not statistically significant and we omitted reporting our results thereof.
We first maximize the likelihood L * 2T to estimate the marginal distributions of financial market returns. As we have shown before, we confirmed that the marginal distributions of market returns have thicker tails than the normal distribution in Table 9. Hence, it may be appropriate to use the GPD in our estimation. It is because the estimates of i (i = 1, 2) are positive and it means that the Frechet-type tail distribution is appropriate The estimated model consists of two dimensions ( d = 2 and p = 3 ), and we take the impact functions c(x) as Case (1) c(x) = 1 and Case (2) c(x) = x . The estimated values of the log-likelihood and AIC are those with the marginal distributions L * 1T . The full likelihood can be calculated using L * 1T and L * 2T . The SD of the estimated coefficients are evaluated numerically using the inverse of the estimated Fisher information matrix.

Model 1
We estimated the intensity function as Again the ML estimates can sometimes be numerically unstable, we set restrictions so that the discounted parameters ij (i, j = 1, 2, 3) have the same value . We show the estimation results in Table 10. Note that in the above table N 1 of Model 1 corresponds to Tokyo and N 2 corresponds to Hong Kong in Tokyo-Hong Kong markets.

Model 2
We estimated the intensity function as  We present our estimation results in Table 11. From our estimated results, we find that Model 2 is better than Model 1 as in Example 1. When comparing Tables 10 and 11, we see several interesting findings. The value of AIC in Model 2 is better than Model 1 as we observed in the Tokyo-New York and Tokyo-London datasets. The estimates of coefficients of past effects are often statistically insignificant in the estimated intensity functions ( 12 and 21 ), while the contemporaneous effects of the co-jump term are statistically significant ( 13 and 23 ). This aspect basically agrees with our motivations for developing the SHPP models.

Non-causality tests
In applying the Granger non-causality test procedure, we set the impact function as c(x) = x . We report our empirical results for the hypothesis H 0 ∶ ij = 0 using LRT statistics.  Japan to Japan US to Japan Japan to Japan UK to Japan For the null-hypothesis H 0 ∶ 13 = 0, LRT statistic based on Tokyo-Hong Kong data was 11.14 and we reject the null-hypothesis. (The upper 95% critical value of 2 (1) is 3.481). Thus, we revealed a significant instantaneous causal relationship between the Japanese financial market and Hong-Kong financial markets.
For testing the null-hypothesis H 0 ∶ 12 = 0, LRT statistics was 0.0, and the null hypothesis was accepted. In addition, for testing the null-hypothesis H 0 ∶ 12 = 0, 13 = 0 , LRT statistic was 11.14; thus the null-hypothesis was rejected.
For the null-hypothesis H 0 ∶ 21 = 0, LRT statistic was 0.006 and we cannot reject the null-hypothesis. (The upper 95% critical value of 2 (1) is 3.481). For testing the null-hypothesis H 0 ∶ 23 = 0, LRT statistic was 2.42, and the null-hypothesis was accepted. Similarly, for the null-hypothesis H 0 ∶ 21 = 0, 23 = 0 LRT statistic was 2.66, and the null-hypothesis could not be rejected.
To summarize our findings in this subsection among Tokyo and Hong Kong financial markets, we found that the simultaneous effects of the two markets are significant, while the effects of past events are rather small.

Further empirical analysis
Here we employ spectral decomposition and RPC as explained in Sect. 3; see Figs. 5, 6 and 7 for the United States, United Kingdom, and Hong Kong, respectively. In the former (two) decompositions, we assume there are no co-jumps, while in the last one co-jumps are permitted. We adopted Models with c ij (x) = x Japan to Japan HK to Japan co−jump to Japan because the resulting models minimize AIC. The graphs are depicted at each frequency, but truncated in x axis because it seems that our empirical data of point processes do not have much information in high frequencies.
In particular, Fig. 5 gives the spectral decomposition based on the estimated intensity model (Model 2) from Japan (Nikkei-225) data, which gives the relative contributions from Japan itself and from US ( S&P 500). Figure 6 gives the spectral decomposition based on the estimated intensity model (Model 2) from Japan (Nikkei-225) data, which gives the relative contributions from Japan itself and from UK (FTSE). Then Fig. 7 gives the spectral decomposition based on the estimated intensity model with co-jumps (Model 2) from Japan (Nikkei-225) data, which gives the relative contributions from Japan itself, Hong Kong (Hansen) and the instantaneous relation.
From these figures, we found that for the relationship between the Tokyo-New York financial markets, self contribution of past events plays a major role, while there is some contribution from New York to Tokyo in the low frequency, which corresponds to the long-run relationship. On the other hand, for the relationship between the Tokyo-Hong Kong financial markets, the instantaneous contribution and the self contribution play major roles in all frequencies. This aspect reflects the fact that we used SHPP models.

Conclusions
In this paper, we developed a new method of econometric analysis of multivariate time series of events and proposed the simultaneous multivariate Hawkes-type point process (SHPP) modeling. Unlike some existing studies, we developed and used new statistical models for simultaneous sudden, large events and delayed events occurring explicitly. Using the SHPP models, we investigated Granger causality and instantaneous Granger causality on several financial markets and economies, and developed bespoke non-causality tests.
By applying GNC and IGNC tests, we revealed the important relationships among major financial markets and several empirical findings. In the Tokyo-New York financial markets, there is a strong unidirectional causation, while in the Tokyo-Hong Kong financial markets the simultaneous effects are dominant.
Several questions remain to be answered. First, although we used Hawkes-type marked point processes, there can be many possible non-linear point processes and Kurisu (2018) discussed one way to justify the use of SHPP models. In economic and financial econometrics, it is standard to handle discrete time series observations in yearly, monthly, weekly, daily, hourly, and per minute terms. Thus, we need a coherent way of investigating abrupt or sudden events, and we propose one way to deal with discrete time series events in this paper. In this respect, it should be interesting to investigate the robustness of our empirical results further. Second, the choice of threshold parameter is an important issue that is related to the relevance of the GPD in SEVT. Since we used a simple threshold parameter, we need a convincing justification on the choice of threshold. Finally, when d > 2 there can be many parameters to be estimated, and the estimated parameters would often be statistically insignificant. This aspect is important when we have the possibility of co-jumps and we used the likelihood and its AIC as a criterion. We need to investigate the problem of choosing statistical point process models further.
These issues are currently under investigation, we shall report our progress in these respects on another occasion. where ( 0 ) is the Fisher information matrix.
Remark A-2 As corollaries of Theorems A.2 and A.3, it is straightforward and standard, but lengthy to give the formal proofs of Theorems 1 and 2 as the non-causality tests we have developed and discussed in Sect. 4.

Remark A-3
As a final remark, we re-emphasize that while Ogata (1978) discussed a set of sufficient conditions for the consistency and asymptotic normality of the ML estimator in one-dimensional self-exciting point processes, we extended his results significantly to the multivariate point processes under a set of weaker conditions. For instance, ( 0 ) is not necessarily a constant matrix and our conditions mean the mixed Gaussian distribution in the present formulation of the Appendix. The limiting 2 property of the statistics is often called the Wilks Property.