{"id":"e86d5f4a-e992-4d14-8e6d-a562c8b1a75d","arxiv_id":"2506.08775","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"A Markovian multivariate Hawkes population process has a closed-form joint transform, with recursive formulas and block-matrix algorithms for all transient and stationary moments.","lead":"This paper derives closed-form formulas for the joint distribution transform and moments of multivariate Hawkes processes with an attached population process, and shows how to compute those moments efficiently. The practical payoff is fast moment-based estimation and covariance analysis for models of contagion in finance, epidemiology, and service systems.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Moment recursions in Eqs. (21)/(23) replace joint moments of jump-size vectors by products of marginal moments, invalidating the claimed generality for dependent intensity jumps.","rationale":"The reader’s weakest assumption was about Markovianity outside the exponential-decay case, which is a scope limitation the paper itself acknowledges. My review found a more direct internal correctness problem: the derivation of the joint moment ODEs (21) and (23) from the transform PDE contains an unjustified factorization of joint moments of the jump-size vector into products of marginal moments. This is load-bearing because the paper’s central claimed contribution is exactly the moment recursions for a general multivariate Markovian Hawkes process, with arbitrary joint distributions of the intensity jump sizes. The error does not affect Theorem 1’s transform characterization, which correctly uses the joint Laplace transform β_j(s) = E[e^{-s^T B_j}], but it invalidates the claimed moment consequences in the stated generality. The numerical experiments all assume independent components of B_j, so the error is not detected there. A concrete analytic check for the bivariate cross-moment ODE shows the coefficient discrepancy explicitly. Because the central claim—valid moment recursions for general jump-size distributions—is false as stated, the appropriate verdict is REJECT, not merely CONDITIONAL. The paper could be repaired by replacing ∏_i E[B_{ij}^{...}] with E[∏_i B_{ij}^{...}] in Eqs. (21) and (23), or by explicitly restricting the model to independent jump-size components, but either change is substantial and would affect the presented results. The reader’s conditional accept was based on secondary issues; the issue identified here is more fundamental and was not flagged in the reader’s verdict.","tokens_in":48336,"tokens_out":18267,"duration_ms":202286,"concrete_test":"Take d=2 with α1=α2=1, λ1=λ2=1, μ1=μ2=1, choose B11 and B21 perfectly correlated with P(B11=B21=2)=P(B11=B21=0)=1/2, and set B12=B22=0.5 deterministic. Verify that Assumption 1 holds. Compute E[λ1(t)λ2(t)] two ways: (a) solve the exact second-moment ODE derived directly from the generator PDE (58) or simulate the SDE; (b) apply the printed Eq. (21). The printed recursion uses E[B11]E[B21]=1 where the correct coefficient is E[B11B21]=2. If the two moment trajectories differ (beyond numerical error), the concern is confirmed.","verdict_should_be":"REJECT","load_bearing_attack":"The central moment recursions are derived from the generator PDE (58) by differentiation, but the passage from Eq. (62) to Eq. (63) in Appendix B.2 factorizes incorrectly. In Eq. (62), the relevant term is E[ z_j ∏_i B_{ij}^{n_{λ,i}-m_i} λ_i^{...} z^{Q_i} ]. Since B_j is independent of (λ,Q), this equals E[∏_i B_{ij}^{n_{λ,i}-m_i}] · E[z_j ∏_i λ_i^{...} z^{Q_i}], so the coefficient is the joint moment E[∏_i B_{ij}^{n_{λ,i}-m_i}]. Eq. (63) instead writes ∏_i E[B_{ij}^{n_{λ,i}-m_i}], which is equal only when the components of B_j are independent across i. For example, with d=2, nλ=(1,1), nQ=0, and j=1, the last sum of Eq. (21) gives E[B11]E[B21] E[λ1], whereas the same term obtained directly from the generator (or from Itô’s formula for d(λ1λ2)) is E[B11B21] E[λ1]. These differ unless B11 and B21 are independent. The paper explicitly allows general joint distributions of the intensity jump sizes (Introduction, Contributions), so Eqs. (21) and (23) are not the claimed general moment recursions. Theorem 1 itself remains correct; the error enters only in the differentiation step used to derive the moment systems. The numerical experiments in Section 7 use independent exponentially distributed marks, which hides the discrepancy.","agreement_with_reader":"disagree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper studies Markovian multivariate Hawkes processes with exponential decay and their induced infinite-server population processes. It derives a closed-form joint transform for (Q(t), λ(t)) via the Markov generator and the method of characteristics (Theorem 1), extends this to two time points (Theorem 2), and then differentiates the transform to obtain linear ODE systems for transient moments (Eq. (21)) and algebraic systems for stationary moments (Eq. (23)). For the bivariate case it develops recursive and nested block-matrix algorithms for moments of arbitrary order, studies the nearly unstable regime under a symmetry assumption (Theorem 3, Corollary 2), and reports numerical experiments comparing the proposed moment computations with finite differences and Monte Carlo simulation.","tokens_in":48661,"tokens_out":5462,"duration_ms":71454,"significance":"If correct, the paper would provide a substantial computational toolkit: exact moment recursions for all transient and stationary cross-moments of Markovian multivariate Hawkes population processes, closed-form autocovariance evaluation via Theorem 2, and an explicit heavy-traffic-type limit for the intensity. The transform results in Theorems 1 and 2 are derived in detail and appear sound. The paper also provides reproducible code, extensive appendices, and careful numerical comparisons. However, the moment recursion claimed for general joint distributions of the intensity jump sizes contains a factorization error that invalidates Eqs. (21) and (23) as stated; this is a load-bearing issue for the paper's central claims, even though the underlying transform theorem remains correct.","major_comments":[{"comment":"The passage from Eq. (62) to Eq. (63) factorizes E[z_j ∏_i B_{ij}^{n_{λ,i}-m_i} λ_i(t)^{...} z^{Q_i(t)}] incorrectly. Since B_j is independent of (λ,Q), the term equals E[∏_i B_{ij}^{n_{λ,i}-m_i}] · E[z_j ∏_i λ_i^{...} z^{Q_i}], so the coefficient is the joint moment E[∏_i B_{ij}^{n_{λ,i}-m_i}], not the product ∏_i E[B_{ij}^{n_{λ,i}-m_i}] written in Eq. (63) and in Eq. (21). These are equal only when the components of B_j are independent across i. The paper explicitly allows general joint distributions of the intensity jump sizes (Eq. (2), Introduction), and Eq. (18) even assumes finiteness of the joint moment E[∏_i B_{ij}^{n_{λ,i}}], so the factorization is not justified. A concrete consequence: for d=2, nλ=(1,1), nQ=0, j=1, the offending term in Eq. (21) gives E[B11]E[B21] E[λ1] instead of the direct-differentiation result E[B11B21] E[λ1]. The same error propagates into the stationary system Eq. (23), the recursive algorithms of Section 4, and the nested block-matrix construction of Section 5. The numerical experiments in Section 7 use independent exponentially distributed marks, which hide the discrepancy. The fix is local in principle—replace the products of marginal B-moments with the joint mixed moments of B_j—but it affects all moment-vector and block-matrix formulas and must be carried through Sections 4, 5, and Appendices C-D.","section":"Appendix B.2, Eqs. (62)-(63) and Eq. (21)"},{"comment":"Because the recursive and nested block-matrix procedures are built directly on Eq. (21), their validity for general dependent marks is currently unsupported. The matrices in Section 5 and Appendix D.4 contain objects such as E[B11]E[B21] and E[B_{11}^2]E[B_{21}] where the correct coefficients are joint moments E[B11B21] and E[B_{11}^2 B_{21}]. The presentation should either state and prove the recursions under an explicit independence assumption on the B_{ij} across i for each j, or replace every occurrence of product marginal B-moments by the corresponding joint moments and update the block-matrix formulas accordingly.","section":"Sections 4-5 and Appendix D"}],"minor_comments":[{"comment":"The parameterization of the limiting Gamma distribution is inconsistent with Theorem 3. From Eq. (49), the Laplace transform is (σ/(σ+s))^{σλ}, which corresponds to a Gamma distribution with shape σλ and rate σ (or scale 1/σ), not scale σ as stated in Corollary 2 and in the sentence defining Γ(r, λ). The covariance formula in Eq. (50) is consistent with the common-component interpretation, but the corollary should also state explicitly that the limiting vector is degenerate in the sense that all components are equal to a single common Gamma random variable.","section":"Corollary 2, Section 6"},{"comment":"In the displayed computation for the sixth term after Eq. (56), the second sum is written with a factor k_j ν_j inside the Laplace transform; this should be ν_j, since the k_j μ_j part is already separated. This is a typographical issue that does not affect the result.","section":"Appendix B.1, sixth term derivation"},{"comment":"The definition of (B_ij(t))_t as 'a sequence of independent random variables' leaves unclear whether independence is assumed across i for fixed j and across j for fixed i. Given that the paper's generality claim depends on the joint law of B_j, it would help to state explicitly what is independent and what may be dependent.","section":"Definition 1, Eq. (2)"}],"recommendation":"major_revision","confidential_remarks":"The core transform results (Theorems 1 and 2) appear correct and are valuable, but the paper's advertised moment recursions for general jump-size distributions contain a genuine factorization error. The numerical experiments, which use independent marks, do not expose the issue. I recommend major revision rather than rejection because the error is localized to the differentiation step and can in principle be repaired by using joint moments of B_j throughout; after such a repair the algorithmic structure may still hold, but the revised manuscript must be re-checked for consistency across Sections 4, 5, and the appendices."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Two things to know. First, Theorem 1 — the closed-form joint transform for the Markovian multivariate Hawkes population process — is genuine, carefully derived, and a real extension of the univariate results. Second, the central moment recursions in Eqs. (21) and (23) are not valid for the general jump-size distributions the paper advertises. In Appendix B.2, the passage from Eq. (62) to Eq. (63) factorizes E[∏_i B_{ij}^{a_i}] as ∏_i E[B_{ij}^{a_i}], which is only true when the components of B_j are independent. The paper never states that assumption, and the Introduction explicitly says general distributions are allowed. The numerical experiments use independent exponential marks, so the error is hidden.\n\nWhat is genuinely good: the transform characterization via the generator PDE and method of characteristics is solid, and the two-time transform in Theorem 2 is a nice addition that enables autocovariance computations. The nested block-matrix construction in Section 5 is a useful computational organization, and the speed/accuracy claims for independent marks are credible from the tables. Theorem 3, the nearly-unstable Gamma limit under symmetry, is a clean and believable extension.\n\nThe soft spots are mostly downstream of the factorization error. The recursive ODE systems, the block-matrix algorithms, and the stationary algebraic systems all inherit the problem. For independent marks they are correct, but the paper needs to say so clearly. Two smaller issues: the relation to the authors' own [17] is left vague, and Corollary 2's Gamma notation is easy to misread (shape/rate versus shape/scale), though the covariance expression is right.\n\nWho should read it: anyone working on Hawkes-fed infinite-server queues who wants exact transforms, and anyone doing moment-based estimation of multivariate Hawkes models — but the latter should wait until the recursion issue is resolved. This paper deserves a serious referee because Theorem 1 alone is worth the referee time. The referee should ask the authors to either state the independence assumption explicitly and check where it enters, or fix the factorization and rerun the numerics. Either way, the paper can be made correct and useful.","headline":"The joint transform is a solid new result, but the moment recursions factorize joint jump-size moments incorrectly and only hold for independent marks.","tokens_in":49190,"tokens_out":3959,"would_cite":true,"duration_ms":47493,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["60G55","60J25","60K25","60E10"],"pacs":[],"model":"deepseek-v4-flash","headline":"For Markovian multivariate Hawkes processes with exponential decay, the joint transform of population and intensity has a closed form, and differentiating it yields exact transient and stationary cross-moments of any order.","keywords":["Hawkes processes","mutual excitation","Markov processes","transform analysis","moment computations","infinite-server queues","population processes","nearly unstable regime"],"falsifier":"Set the bivariate parameters as in Section 7.2, solve the ODE system (9) for a nontrivial pair $(s,z)$, and compare the value of $\\zeta_0(t,s,z)$ from Eq. (8) with a Monte Carlo estimate of $\\mathbb{E}[z_1^{Q_1(t)}z_2^{Q_2(t)}e^{-s_1\\lambda_1(t)-s_2\\lambda_2(t)}]$ using enough replications that the Monte Carlo standard error is negligible; a persistent gap beyond that error would refute the closed-form transform, and agreement would confirm it.","tokens_in":48113,"feed_emoji":"📈","tokens_out":8596,"duration_ms":94695,"temperature":0.7,"pith_summary":"This paper proves that for Markovian multivariate Hawkes processes—point processes whose intensities decay exponentially and jump at every event—the joint transform of the population size and the current intensity is available in closed form as a product of exponentials whose arguments solve a small system of ordinary differential equations. Because the transform is bijective to the joint distribution, differentiating it yields exact analytic expressions for transient and stationary cross-moments of any order, together with auto- and cross-covariance functions, rather than simulation-based estimates. The paper organizes the moment relations into a nested sequence of block lower triangular matrices, which makes the computation fast and, for transient moments, insensitive to the time horizon. It also shows that in a symmetric parameterization the stationary intensity, scaled by the distance to instability, converges to a multivariate Gamma law. If these results are correct, moment-based inference and comparative statics for contagious-event models become practical in settings where finite differences and Monte Carlo simulation were the only options.","feed_headline":"One transform yields all moments of Hawkes populations","feed_subtitle":"Exact transient and stationary moment formulas replace simulation for mutually-exciting event models.","key_machinery":"The load-bearing object is the joint transform $\\zeta_{t_0}(t,s,z)$, a combined z-transform in $Q(t)$ and Laplace transform in $\\lambda(t)$ that is bijective to the joint law of the population and intensity processes. The derivation writes the Markov generator of $(Q,\\lambda)$ as a first-order PDE, (58), and solves it by the method of characteristics; the characteristic equations are precisely the ODE system (9). The moment machinery is repeated differentiation of this transform, yielding the moment recursions (21) and (23). The computational shortcut is the nested sequence of block lower triangular matrices $A_n$ and $F_n$ built from the blocks $M^{(k,n-k)}$, whose invertibility turns the moment ODEs into closed-form matrix-exponential expressions (40) and (44).","core_discovery":"The central claim is Theorem 1: for a $d$-component Hawkes process with exponential decay rates $\\alpha_i$ and independent jump marks $B_{ij}$, and for the induced infinite-server population $Q(t)$, the conditional joint transform $\\zeta_{t_0}(t,s,z)=\\mathbb{E}[\\prod_{i=1}^d z_i^{Q_i(t)}e^{-s_i\\lambda_i(t)}]$ is given in closed form by (8), where the functions $b_{z_j}$ and $\\tilde{s}_j$ solve the explicit ODE system (9). Repeated differentiation of this transform produces the linear ODE recursion (21) for transient reduced moments and the algebraic recursion (23) for stationary moments, both of total order at most the target order, so all cross-moments can be solved recursively. In the bivariate setting the recursions are reorganized as nested block lower triangular matrices, and Proposition 4 gives transient moments up to order $n$ as a single matrix-exponential formula whose run time does not grow with the horizon. Separately, Theorem 3 establishes that, under a symmetric parameterization and with $\\theta = \\frac{1}{\\alpha}\\sum_i\\mathbb{E}[B_i]\\uparrow 1$, the Laplace transform of the scaled stationary intensity obeys $\\lim_{\\theta\\uparrow1}\\mathcal{T}\\{\\lambda\\}(s(1-\\theta))=(\\sigma/(\\sigma+s))^{\\sigma\\lambda}$, a multivariate Gamma limit.","pith_inferences":["An implication the paper leaves implicit is that the transform-differentiation recipe depends only on the generator being a polynomial in the state variables, so the same scheme should extend to other Markovian intensity models, such as Markov-modulated or affine jump intensities, without reworking the recursion.","For non-exponential decay kernels the Markov property fails and Theorem 1 should not be applied; a natural test is to approximate a general kernel by a phase-type (sum-of-exponentials) kernel and compare the moment formulas against a cluster-representation simulation of the true Hawkes process.","The nearly-unstable Gamma limit suggests a further, unproved consequence: in the same parameterization the scaled population process $Q(t)$ may also possess a tractable heavy-traffic limit, by analogy with infinite-server queues in overload.","Because Theorem 2 gives the full two-time transform, spectral or frequency-domain inference for Hawkes populations is a plausible next step, though the paper does not explore it."],"forward_implications":["Transient and stationary cross-moments of any order, as well as auto- and cross-covariances, can be evaluated exactly by solving linear ODEs or linear algebraic systems, without simulation error.","The nested block-matrix construction makes the cost of transient moment evaluation essentially independent of the time horizon, because only a matrix exponential of fixed dimension is involved.","Moment-based estimation and comparative statics for multivariate Hawkes models become practical, since each evaluation of a collection of moments is near-instant and exact.","In the nearly unstable symmetric regime, the stationary intensity, scaled by $1-\\theta$, converges to a Gamma law with covariance $\\lambda/\\sigma$ between every pair of components.","The two-time transform of Theorem 2 supplies the processes' autocovariance functions directly, enabling likelihood-free inference that uses temporal dependence information."],"supporting_citations":[{"why":"Supplies the Markov property for exponentially decaying self-exciting processes, which is the premise for the transform PDE.","marker":"[21]"},{"why":"Establishes the univariate infinite-server queue with Hawkes input whose moment recursions this paper generalizes to the multivariate setting.","marker":"[18]"},{"why":"Provides the matrix-calculus perspective for moments of Markov processes that the nested block-matrix construction adapts.","marker":"[11]"},{"why":"Introduces mutually exciting point processes and their spectra, the model class under study.","marker":"[14]"},{"why":"Gives the standard conditional-intensity characterization of point processes used in Definition 1.","marker":"[8]"},{"why":"Supplies the stability condition for nonlinear Hawkes processes that underlies Assumption 1.","marker":"[3]"},{"why":"Gives the cluster representation of general multivariate Hawkes populations, used as the benchmark extension and basis for simulation.","marker":"[17]"}],"fun_headline_variants":["Exact moments for multivariate Hawkes processes","One formula yields every Hawkes moment","Closed-form moments for Hawkes population models","Efficient moment evaluation for Hawkes processes","All Hawkes moments from a single transform"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The whole construction rests on the joint process $(N(t),\\lambda(t))$ being Markov, which holds only when the excitation decay functions are exponential and the jump marks are independent and identically distributed; drop that and the closed form in Theorem 1 is no longer claimed.","fun_headline_variants_meta":{"raw":{"variants":["Exact moments for multivariate Hawkes processes","One formula yields every Hawkes moment","Closed-form moments for Hawkes population models","Efficient moment evaluation for Hawkes processes","All Hawkes moments from a single transform"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000242,"raw_usage":{"total_tokens":1531,"prompt_tokens":959,"completion_tokens":572,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":575,"completion_tokens_details":{"reasoning_tokens":507}},"tokens_in":575,"tokens_out":572,"duration_ms":6633,"temperature":1.0,"reasoning_tokens":507,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-07T05:02:55.156078+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Set the bivariate parameters as in Section 7.2, solve the ODE system (9) for a nontrivial pair $(s,z)$, and compare the value of $\\zeta_0(t,s,z)$ from Eq. (8) with a Monte Carlo estimate of $\\mathbb{E}[z_1^{Q_1(t)}z_2^{Q_2(t)}e^{-s_1\\lambda_1(t)-s_2\\lambda_2(t)}]$ using enough replications that the Monte Carlo standard error is negligible; a persistent gap beyond that error would refute the closed-form transform, and agreement would confirm it.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the Markov property for exponentially decaying self-exciting processes, which is the premise for the transform PDE."},{"cited_title":"T., Saxena, M., Boxma, O","cited_arxiv_id":null,"evidence_quote":"Establishes the univariate infinite-server queue with Hawkes input whose moment recursions this paper generalizes to the multivariate setting."},{"cited_title":"& Pender, J","cited_arxiv_id":null,"evidence_quote":"Provides the matrix-calculus perspective for moments of Markov processes that the nested block-matrix construction adapts."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Introduces mutually exciting point processes and their spectra, the model class under study."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Gives the standard conditional-intensity characterization of point processes used in Definition 1."},{"cited_title":"& Massouli´e, L","cited_arxiv_id":null,"evidence_quote":"Supplies the stability condition for nonlinear Hawkes processes that underlies Assumption 1."}],"review_version":1}