3.1. The moments of general polynomial Markov processes

The connection between matryoshkan matrices and Markov processes is built upon a key tool for Markov processes, the infinitesimal generator. For a detailed introduction to infinitesimal generators and their use in studying Markov processes, see e.g. [Reference Ethier and Kurtz23]. For a Markov process $X_t$ on a state space $\mathbb{S}$ , the infinitesimal generator on a function $f \,:\, \mathbb{S} \to \mathbb{R}$ is defined as

\begin{equation*}\mathcal{L}\,f(x)=\lim_{\tau \to 0}\frac{{\mathbb{E}\!\left[f(X_\tau) \mid X_0 = x\right]} - f(x)}{\tau}.\end{equation*}

In our context and in many others, the power of the infinitesimal generator comes through use of Dynkin’s formula, which gives us that

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[f(X_t)\right]}={\mathbb{E}\!\left[\mathcal{L}\,f(X_t)\right]}.\end{equation*}

To study the moments of a Markov process, we are interested in functions f that are polynomials. Let us suppose now that $\mathcal{L}x^n$ for any $n \in \mathbb{Z}^+$ is polynomial in the lower powers of x for a given Markov process $X_t$ , meaning that $X_t$ is a one-dimensional polynomial process. Then we can write

\begin{equation*}\mathcal{L}X_t^n=\theta_{0,n}+\sum_{i=1}^n \theta_{i,n} X_t^i,\end{equation*}

which implies that the differential equation for the nth moment of this process is

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t} {\mathbb{E}\!\left[X_t^n\right]}=\theta_{0,n} + \sum_{i=1}^n \theta_{i,n} {\mathbb{E}\big[X_t^i\big]},\end{equation*}

for some collection of constants $\theta_{0,n}$ , $\theta_{1,n}$ , $\theta_{n,n}$ . Thus, to solve for the nth moment of $X_t$ we must first solve for all the moments of lower order. We can also observe that to solve for the $(n-1)$ th moment we must have all the moments below it. In comparing these systems of differential equations, we can see that all of the equations in the system for the $(n-1)$ th moment are also in the system for the nth moment. No coefficients are changed in any of these lower moment equations; the only difference between the two systems is the inclusion of the differential equation for the nth moment in its own system. Hence, the nesting matryoshkan structure arises. By expressing each system of linear ODEs in terms of a vector of moments, a matrix of coefficients, and a vector of shift terms, we can use these matrix sequences to capture how one moment’s system encapsulates all the systems below it. This observation then allows us to calculate all the moments of the process in closed form, as we now show in Theorem 1. In steady state, this also implies that the $(n+1)$ th moment can be computed as a linear combination of moments 1 through n.

Theorem 1. Let $X_t$ be a Markov process such that the time derivative of its nth moment can be written as

(7) \begin{align}\frac{\textrm{d}}{\textrm{d}t} {\mathbb{E}\!\left[X_t^n\right]}&=\theta_{0,n} + \sum_{i=1}^n \theta_{i,n} {\mathbb{E}\big[X_t^i\big]},\end{align}

for any $n \in \mathbb{Z}^+$ , where $t \geq 0$ and $\theta_{i, n} \in \mathbb{R}$ for all $i \leq n$ . Let $\boldsymbol{\Theta}_n \in \mathbb{R}^{n \times n}$ be defined recursively by

(8) \begin{align}{\boldsymbol{\Theta}}_n=\left[\begin{array}{c@{\quad}c}\boldsymbol{\Theta}_{n-1} & \textbf{0}_{n-1 \times 1} \\[5pt]\boldsymbol{\theta}_n & \theta_{n,n}\end{array}\right],\end{align}

where $\boldsymbol{\theta}_n = [\theta_{1, n}, \, \dots, \, \theta_{n-1, n}]$ and $\boldsymbol{\Theta}_1 = \theta_{1,1}$ . Furthermore, let $\boldsymbol{\theta}_{0,n} = [\theta_{0,1}, \, \dots, \, \theta_{0,n}]^{\textrm{T}}$ . Then, if $\theta_{k,k} \ne 0$ for all $k \in \{1, \dots , n\}$ are distinct, the nth moment of $X_t$ can be expressed as

(9) \begin{align}{\mathbb{E}\!\left[X_t^n\right]}&=\boldsymbol{\theta}_{n}\!\left(\boldsymbol{\Theta}_{n-1} - \theta_{n,n}\textbf{I}\right)^{-1}\big(e^{\boldsymbol{\Theta}_{n-1} t} - e^{\theta_{n,n}t}\textbf{I}\big)\left(\textbf{x}_{n-1}(x_0)+\frac{\boldsymbol{\theta}_{0,n-1}}{\theta_{n,n}}\right)\nonumber\\&\quad+x_0^n e^{\theta_{n,n}t}-\frac{\theta_{0,n}\big(1-e^{\theta_{n,n}t}\big)}{\theta_{n,n}}+\frac{\boldsymbol{\theta}_{n}}{\theta_{n,n}}\boldsymbol{\Theta}_{n-1}^{-1}\big(\textbf{I}-e^{\boldsymbol{\Theta}_{n-1}t}\big)\boldsymbol{\theta}_{0,n-1},\end{align}

where $x_0$ is the initial value of $X_t$ and where $\textbf{x}_n(a) \in \mathbb{R}^n$ is such that $\left(\textbf{x}_n(a)\right)_i = a^i$ . If $X_t$ has a stationary distribution and the nth steady-state moment ${\mathbb{E}\!\left[X_\infty^n\right]}$ exists, it is given by

(10) \begin{align}{\mathbb{E}\!\left[X_\infty^n\right]}&=\frac{1}{\theta_{n,n}}\big(\boldsymbol{\theta}_{n}\boldsymbol{\Theta}_{n-1}^{-1}\boldsymbol{\theta}_{0,n-1}-\theta_{0,n}\big),\end{align}

and these steady-state moments satisfy the recursive relationship

(11) \begin{align}{\mathbb{E}\big[X_\infty^{n+1}\big]}&=-\frac{1}{\theta_{n+1,n+1}}\!\left(\boldsymbol{\theta}_{n+1}\textbf{s}_{n}^X(\infty)+\theta_{0,n+1}\right),\end{align}

where $\textbf{s}_n^X(\infty) \in \mathbb{R}^n$ is the vector of steady-state moments such that $\left(\textbf{s}_n^X(\infty)\right)_i = {\mathbb{E}\!\left[X_\infty^i\right]}$ .

Proof. Using the definition of $\boldsymbol{\Theta}_n$ in Equation (8), Equation (7) gives rise to the system of ODEs given by

(12) \begin{align}\frac{\textrm{d}}{\textrm{d}t}\textbf{s}_n^X(t)=\boldsymbol{\Theta}_n \textbf{s}_n^X(t)+\boldsymbol{\theta}_{0,n},\end{align}

where $\textbf{s}_n^X(t) \in \mathbb{R}^n$ is the vector of transient moments at time $t \geq 0$ such that $\left(\textbf{s}_n^X(t)\right)_i = {\mathbb{E}\big[X_t^i\big]}$ for all $1 \leq i \leq n$ . We can observe that by definition the matrices $\boldsymbol{\Theta}_n$ form a matryoshkan sequence, and thus by Lemma 1 we achieve the stated transient solution. To prove the steady-state solution, we can first note that if the process has a steady-state distribution then the vector $\textbf{s}_n^X(\infty) \in \mathbb{R}^n$ defined by $\left(\textbf{s}_n^X(\infty)\right)_i = {\mathbb{E}\!\left[X_\infty^i\right]}$ will satisfy

(13) \begin{align}0=\boldsymbol{\Theta}_n \textbf{s}_n^X(\infty)+\boldsymbol{\theta}_{0,n},\end{align}

as this is the equilibrium solution to the differential equation corresponding to each of the moments. This system has a unique solution, since $\boldsymbol{\Theta}_n$ is nonsingular owing to the assumption that the diagonal values are unique and nonzero. Using Proposition 1, we find the nth moment by

\begin{equation*}{\mathbb{E}\!\left[X_\infty^n\right]}=-\textbf{v}_n^{\textrm{T}} \boldsymbol{\Theta}_n^{-1} \boldsymbol{\theta}_{0,n}=\left[\begin{array}{c@{\quad}c}\frac{1}{\theta_{n,n}}\boldsymbol{\theta}_{n}\boldsymbol{\Theta}_{n-1}^{-1} & -\frac{1}{\theta_{n,n}}\end{array}\right]\left[\begin{array}{c}\boldsymbol{\theta}_{0,n-1} \\[5pt]\theta_{0,n}\end{array}\right]=\frac{1}{\theta_{n,n}}\Big(\boldsymbol{\theta}_{n}\boldsymbol{\Theta}_{n-1}^{-1}\boldsymbol{\theta}_{0,n-1}-\theta_{0,n}\Big),\end{equation*}

which completes the proof of Equation (10). To conclude, one can note that each line of the linear system in Equation (13) implies the stated recursion.

3.2. Application to Hawkes process intensities

For our first example of this method let us turn to our motivating application, the Markovian Hawkes process intensity. Via [Reference Hawkes33], this process is defined as follows. Let $\lambda_t$ be a stochastic arrival process intensity such that

\begin{equation*}\lambda_t=\lambda^* + (\lambda_0 - \lambda^*) e^{-\beta t} + \int_0^t \alpha e^{-\beta (t-s)} \textrm{d}N_s=\lambda^* + (\lambda_0 - \lambda^*) e^{-\beta t} + \sum_{i=1}^{N_t} \alpha e^{-\beta (t - A_i)},\end{equation*}

where $\big\{A_i \mid i \in \mathbb{Z}^+\big\}$ is the sequence of arrival epochs in the point process $N_t$ , with

\begin{equation*}\mathbb{P}\!\left( {N_{t + s} - N_t = 0 \mid \mathcal{F}_t} \right)=\mathbb{P}\!\left( {N_{t + s} - N_t = 0 \mid \lambda_t} \right)=e^{-\int_0^s \lambda_{t + u} \textrm{d}u},\end{equation*}

where $\mathcal{F}_t$ is the filtration generated by the history of $\lambda_t$ up to time t. Let us note that the Hawkes process as defined in [Reference Hawkes33] is not in general a Markov process. If we restrict the excitation kernel to being exponential, then we obtain a Markovian process. In what remains, we will assume that the excitation kernel is exponential. Moreover, we will assume that $\beta > \alpha > 0$ , so that the process has a stationary distribution, and we will also let $\lambda^* > 0$ and $\lambda_0 > 0$ . Note that the process behaves as follows: at arrivals $\lambda_t$ increases by $\alpha$ , and in the interims it decays exponentially at rate $\beta$ towards the baseline level $\lambda^*$ . Thus, $(\lambda_t, N_t)$ is referred to as a self-exciting point process, as the occurrence of an arrival increases the intensity and thus increases the likelihood that another arrival will occur soon afterwards. Because the intensity $\lambda_t$ forms a Markov process, we can follow the literature (see, e.g., Equation (7) of [Reference Da Fonseca and Zaatour15]) and write its infinitesimal generator for a (sufficiently regular) function $f\,:\,\mathbb{R}^+ \to \mathbb{R}$ as follows:

\begin{equation*}\mathcal{L}\,f(\lambda_t) = \lambda_t \!\left(\,f(\lambda_t + \alpha) - f(\lambda_t)\right) - \beta \!\left(\lambda_t - \lambda^*\right) f^{\prime}(\lambda_t),\end{equation*}

where f ^′(x) is the first derivative of $f({\cdot})$ evaluated at x. Note that this expression showcases the process dynamics that we have described, as the first term on the right-hand side corresponds to the product of the arrival rate and the change in the process when an arrival occurs, while the second term captures the decay.

To obtain the nth moment we must consider $f({\cdot})$ of the form $f(x) = x^n$ . In the simplest case, when $n = 1$ , this formula yields an ODE for the mean, which can be written as

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[\lambda_t\right]}=\alpha{\mathbb{E}\!\left[\lambda_t\right]} - \beta\!\left({\mathbb{E}\!\left[\lambda_t\right]} - \lambda^*\right)=\beta \lambda^*-(\beta - \alpha) {\mathbb{E}\!\left[\lambda_t\right]}.\end{equation*}

By comparison, for the second moment at $n = 2$ we have

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\big[\lambda_t^2\big]}={\mathbb{E}\!\left[\lambda_t \big((\lambda_t + \alpha)^2 - \lambda_t^2\big) - 2 \beta \lambda_t \big(\lambda_t - \lambda^*\big)\right]}=\big(2 \beta \lambda^* + \alpha^2\big){\mathbb{E}\!\left[ \lambda_t \right]}-2(\beta-\alpha){\mathbb{E}\big[\lambda_t^2\big]},\end{equation*}

and we can note that while the ODE for the mean is autonomous, the second moment equation depends on both the mean and the second moment. Thus, to solve for the second moment we must also solve for the mean, leading us to the following system of differential equations:

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}\!\left[\begin{array}{c}{\mathbb{E}\!\left[\lambda_t\right]}\\[5pt] {\mathbb{E}\!\left[\lambda_t^2\right]}\\[2pt]\end{array}\right]=\left[\begin{array}{c@{\quad}c}-(\beta - \alpha) & 0 \\[5pt]2\beta\lambda^* + \alpha^2 & -2(\beta-\alpha)\end{array}\right]\left[\begin{array}{c}{\mathbb{E}\!\left[\lambda_t\right]}\\[5pt]{\mathbb{E}\!\left[\lambda_t^2\right]}\\[2pt]\end{array}\right]+\left[\begin{array}{c}\beta\lambda^* \\[3pt] 0\end{array}\right].\end{equation*}

Moving on to the third moment, the infinitesimal generator formula yields

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\big[\lambda_t^3\big]}&={\mathbb{E}\!\left[\lambda_t\big((\lambda_t + \alpha)^3 - \lambda_t^3\big) - 3\beta \lambda_t^{2}\!\left(\lambda_t - \lambda^*\right)\right]}\\&=\alpha^{3}{\mathbb{E}\!\left[\lambda_t \right]}+3\big(\beta\lambda^* + \alpha^2\big){\mathbb{E}\big[\lambda_t^{2}\big]}-3(\beta - \alpha) {\mathbb{E}\big[\lambda_t^{3} \big]},\end{align*}

and hence we see that this ODE now depends on all of the first three moments. Thus, to solve for ${\mathbb{E}\!\left[\lambda_t^3\right]}$ we need to solve the system of ODEs

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}\!\left[\begin{array}{c}{\mathbb{E}\!\left[\lambda_t\right]}\\[5pt]{\mathbb{E}\!\left[\lambda_t^2\right]}\\[5pt] {\mathbb{E}\!\left[\lambda_t^3\right]}\\[2pt]\end{array}\right]=\left[\begin{array}{c@{\quad}c@{\quad}c}-(\beta - \alpha) & 0 & 0\\[5pt]2\beta\lambda^* + \alpha^2 & -2(\beta-\alpha) & 0 \\[5pt]\alpha^3 & 3\big(\beta\lambda^* + \alpha^2\big) & - 3(\beta-\alpha)\end{array}\right]\left[\begin{array}{c}{\mathbb{E}\!\left[\lambda_t\right]}\\[5pt] {\mathbb{E}\!\left[\lambda_t^2\right]}\\[5pt] {\mathbb{E}\!\left[\lambda_t^3\right]}\\[2pt]\end{array}\right]+\left[\begin{array}{c}\beta\lambda^* \\[5pt] 0 \\[5pt] 0\end{array}\right],\end{equation*}

and this now suggests the matryoshkan structure of these process moments: we can note that the system for the second moment is nested within the system for the third moment. That is, the matrix for the three-dimensional system contains the two-dimensional system in its upper left-hand block, just as the vector of the first three moments has the first two moments in its first two coordinates. In general, we can see that the nth moment will satisfy the ODE given by

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[\lambda_t^n\right]}&={\mathbb{E}\big[\lambda_t\!\left((\lambda_t + \alpha)^n - \lambda_t^n\right) - n\beta \lambda_t^{n-1}\!\left(\lambda_t - \lambda^*\right)\!\big]}\\&=\sum_{k=1}^{n}\!\left(\begin{array}{c}n\\k - 1\end{array}\right)\alpha^{n-k+1}{\mathbb{E}\big[\lambda_t^{k} \big]}-n\beta {\mathbb{E}\!\left[\lambda_t^{n} \right]}+n\beta\lambda^* {\mathbb{E}\big[\lambda_t^{n-1}\big]},\end{align*}

where we have simplified by use of the binomial theorem. Thus, the system of differential equations needed to solve for the nth moment uses the matrix from the $(n-1)$ th system augmented below by the row

(14) \begin{align}\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}\alpha^n & n\alpha^{n-1},& \scriptsize\left(\begin{array}{c}n\\2\\\end{array}\right)\alpha^{n-2},& \dots,& \scriptsize\left(\begin{array}{c}n\\n - 3\\\end{array}\right)\alpha^3,& \scriptsize\left(\begin{array}{c}n\\n - 2\\\end{array}\right)\alpha^2 + n \beta\lambda^*,& -n(\beta-\alpha)\end{array}\right],\end{align}

and buffered on the right by a column of zeros. To collect these coefficients into a coherent structure, let us define the matrix $\boldsymbol{\mathcal{P}}_n(a) \in \mathbb{R}^{n \times n}$ for $a \in \mathbb{R}$ such that

(15) \begin{align}\left(\boldsymbol{\mathcal{P}}_n(a)\right)_{i,j}=\begin{cases}\scriptsize\big(\begin{array}{c}i\\[-3pt]\,j - 1\end{array}\big) a^{i-j+1}, & i \geq j, \\[5pt]0, & i < j .\end{cases}\end{align}

If we momentarily disregard the terms with $\beta$ in the general augment row in Equation (14), one can observe that the remaining terms in this vector are given by the bottom row of the matrix $\boldsymbol{\mathcal{P}}_n(\alpha)$ . Furthermore, by definition $\big\{\boldsymbol{\mathcal{P}}_n(a) \mid n \in \mathbb{Z}^+\big\}$ forms a matryoshkan matrix sequence. We can also note that $\boldsymbol{\mathcal{P}}_n(a)$ can be equivalently defined as

\begin{equation*}\boldsymbol{\mathcal{P}}_n(a)=\sum_{k=1}^na\!\left[\begin{array}{c@{\quad}c}\textbf{0}_{n-k \times n-k} & \textbf{0}_{n-k \times k}\\[6pt]\textbf{0}_{k \times n-k} & \textbf{L}_k(a)\end{array}\right],\end{equation*}

where $\textbf{L}_k(a) = e^{a\textbf{diag}\!\left(1:k-1,-1\right)}$ is the kth lower triangular Pascal matrix; i.e. the nonzero terms in $\textbf{L}_k(1)$ yield the first k rows of Pascal’s triangle. Alternatively, $\boldsymbol{\mathcal{P}}_n(a)$ can be found by creating a lower triangular matrix from the strictly lower triangular values in $\textbf{L}_{n+1}(a)$ . For these reasons, we refer to the sequence of $\boldsymbol{\mathcal{P}}_n(a)$ as matryoshkan Pascal matrices. For brief overviews and beautiful properties of Pascal matrices, see [Reference Brawer and Pirovino8, Reference Call and Velleman9, Reference Zhang49, Reference Edelman and Strang21]. As we have seen in the preceding derivation, matryoshkan Pascal matrices arise naturally in using the infinitesimal generator to calculate moments of Markov processes. This follows from the application of the binomial theorem to jump terms. Now, in the case of the Markovian Hawkes process intensity we find closed-form expressions for all transient moments in Corollary 1.

Corollary 1. Let $\lambda_t$ be the intensity of a Hawkes process with baseline intensity $\lambda^* > 0$ , intensity jump $\alpha > 0$ , and decay rate $\beta > \alpha$ . Then the nth moment of $\lambda_t$ is given by

(16) \begin{align}{\mathbb{E}\!\left[\lambda_t^n\right]}&=\textbf{m}^\lambda_{n}\!\left(\textbf{M}^\lambda_{n-1} + n(\beta-\alpha)\textbf{I}\right)^{-1}\big(e^{\textbf{M}^\lambda_{n-1} t} - e^{-n(\beta-\alpha)t}\textbf{I}\big)\left(\textbf{x}_{n-1}(\lambda_0)-\frac{\beta \lambda^* \textbf{v}_1 }{n(\beta-\alpha)}\right)\nonumber\\&\quad+\lambda_0^n e^{-n(\beta-\alpha)t}+\mathbb{I}_{\{n = 1\}}\frac{\beta\lambda^*}{\beta - \alpha}\big(1 - e^{-(\beta - \alpha) t}\big)\nonumber\\&\quad-\frac{\beta\lambda^*}{n(\beta-\alpha)}\textbf{m}^\lambda_{n}\!\left(\textbf{M}^\lambda_{n-1}\right)^{-1}\Big(\textbf{I}-e^{\textbf{M}^\lambda_{n-1} t}\Big)\textbf{v}_{1},\end{align}

for all $t \geq 0$ and $n \in \mathbb{Z}^+$ , where $\textbf{v}_1 \in \mathbb{R}^n$ is the unit vector in the first coordinate, $\textbf{M}_n^\lambda = \beta \lambda^* \textbf{diag}\!\left(2\,:\,n,-1\right) - \beta \textbf{diag}\!\left(1\,:\,n\right) + \boldsymbol{\mathcal{P}}_n(\alpha)$ , $\textbf{m}_n^{\lambda}=\left[\left(\textbf{M}_n^\lambda\right)_{n,1},\,\dots,\,\left(\textbf{M}_n^\lambda\right)_{n,n-1}\right]$ is given by

\begin{align*}\left(\textbf{m}_n^{\lambda}\right)_{j}&=\begin{cases}\scriptsize\big(\begin{array}{c}n\\j - 1\end{array}\big)\alpha^{n-j+1}& if\ \ j < n-1,\\[8pt]\scriptsize\big(\begin{array}{c}n\\n - 2\end{array}\big) \alpha^{2} + n\beta \lambda^*& if\ \ j =n-1,\end{cases}\end{align*}

and $\textbf{x}_n(a) \in \mathbb{R}^n$ is such that $\left(\textbf{x}_n(a)\right)_i = a^i$ . In steady state, the nth moment of $\lambda_t$ is given by

(17) \begin{align}\lim_{t \to \infty}{\mathbb{E}\!\left[\lambda_t^n\right]}=-\frac{\beta\lambda}{n(\beta - \alpha)}\textbf{m}_n^\lambda \big(\textbf{M}_{n-1}^\lambda\big)^{-1}\textbf{v}_1,\end{align}

for $n \geq 2$ , with $\lim_{t \to \infty} {\mathbb{E}\!\left[\lambda_t\right]} = \frac{\beta \lambda^*}{\beta - \alpha}$ . Moreover, the $(n+1)$ th steady-state moment of the Hawkes process intensity is given by the recursion

(18) \begin{align}\lim_{t \to \infty}{\mathbb{E}\big[\lambda_t^{n+1}\big]} = \frac{1}{(n+1)(\beta-\alpha)}\textbf{m}_{n+1}^\lambda \textbf{s}_n^\lambda ,\end{align}

for all $n \in \mathbb{Z}^+$ , where $\textbf{s}_n^\lambda \in \mathbb{R}^n$ is the vector of steady-state moments defined so that $\left(\textbf{s}_n^\lambda\right)_i = \lim_{t \to \infty}{\mathbb{E}\!\left[\lambda_t^i\right]}$ for $1 \leq i \leq n$ .

3.3. Application to shot noise processes

As a second example of calculating moments through matryoshkan matrices, consider a Markovian shot noise process; see e.g. [Reference Daley and Vere-Jones16, Chapter 6] for an introduction. That is, let $\psi_t$ be defined so that

\begin{equation*}\psi_t=\sum_{i=1}^{N_t} J_i e^{-\beta (t - A_i)},\end{equation*}

where $\beta > 0$ , $\big\{J_i \mid i \in \mathbb{Z}^+\big\}$ is a sequence of independent and identically distributed positive random variables with cumulative distribution function $F_J({\cdot})$ , $N_t$ is a Poisson process at rate $\lambda > 0$ , and $\big\{A_i \mid i \in \mathbb{Z}^+\big\}$ is the sequence of arrival times in the Poisson process. These dynamics yield the following infinitesimal generator:

\begin{equation*}\mathcal{L}\,f(\psi_t)=\lambda \int_0^\infty \!\left(\,f(\psi_t + j) - f(\psi_t)\right) \textrm{d}F_J(\,j) - \beta \psi_t f^{\prime}(\psi_t),\end{equation*}

which can be simplified directly from Equation (2.2) of [Reference Dassios and Zhao17]. We can note that this is similar to the Hawkes process discussed in Subsection 3.2, as the right-hand side contains a term for jumps and a term for exponential decay. However, this infinitesimal generator formula also shows key differences between the two processes, as the jumps in the shot noise process are of random size and they occur at the fixed, exogenous rate $\lambda > 0$ . Supposing the mean jump size is finite, this now yields that the mean satisfies the ODE

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[\psi_t\right]}=\lambda {\mathbb{E}\!\left[J_1\right]} - \beta {\mathbb{E}\!\left[\psi_t\right]},\end{equation*}

whereas if ${\mathbb{E}\big[J_1^2\big]} < \infty$ , the second moment of the shot noise process is given by the solution to

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\big[\psi_t^2\big]}={\mathbb{E}\!\left[\lambda \big( (\psi_t + J_1)^2 - \psi_t^2 \big) - 2 \beta \psi_t^2\right]}=\lambda {\mathbb{E}\big[J_1^2\big]} + 2\lambda {\mathbb{E}\!\left[J_1\right]} {\mathbb{E}\!\left[\psi_t\right]} - 2 \beta {\mathbb{E}\big[\psi_t^2\big]},\end{equation*}

which depends on both the second moment and the mean. This gives rise to the linear system of differential equations

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}\!\left[\begin{array}{c}{\mathbb{E}\!\left[\psi_t\right]}\\[5pt] {\mathbb{E}\big[\psi_t^2\big]}\\[2pt] \end{array}\right]=\left[\begin{array}{c@{\quad}c}-\beta & 0 \\[5pt] 2\lambda{\mathbb{E}\!\left[J_1\right]} & -2\beta\end{array}\right]\left[\begin{array}{c}{\mathbb{E}\!\left[\psi_t\right]}\\[5pt] {\mathbb{E}\big[\psi_t^2\big]}\\[2pt]\end{array}\right]+\left[\begin{array}{c}\lambda {\mathbb{E}\!\left[J_1\right]} \\[5pt] \lambda {\mathbb{E}\big[J_1^2\big]}\\[2pt]\end{array}\right],\end{equation*}

and by observing that the differential equation for the third moment depends on the first three moments if the third moment of the jump size is finite,

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\big[\psi_t^3\big]} =\, & {\mathbb{E}\!\left[\lambda \!\left( (\psi_t + J_1)^3 - \psi_t^3 \right) - 3 \beta \psi_t^3\right]}=\lambda {\mathbb{E}\big[J_1^3\big]} + 3\lambda {\mathbb{E}\big[J_1^2\big]} {\mathbb{E}\!\left[\psi_t\right]}\\& + 3\lambda{\mathbb{E}\!\left[J_1\right]} {\mathbb{E}\big[\psi_t^2\big]} - 3 \beta {\mathbb{E}\big[\psi_t^3\big]},\end{align*}

we can see that the system for the first two moments is again contained in the system for the first three moments:

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}\!\left[\begin{array}{c}{\mathbb{E}\!\left[\psi_t\right]}\\[5pt] {\mathbb{E}\big[\psi_t^2\big]}\\[5pt] {\mathbb{E}\big[\psi_t^3\big]}\end{array}\right]=\left[\begin{array}{c@{\quad}c@{\quad}c}-\beta & 0 & 0 \\[5pt]2\lambda{\mathbb{E}\!\left[J_1\right]} & -2\beta & 0 \\[5pt]3\lambda{\mathbb{E}\big[J_1^2\big]} & 3\lambda{\mathbb{E}\!\left[J_1\right]} & - 3\beta\\[3pt]\end{array}\right]\left[\begin{array}{c}{\mathbb{E}\!\left[\psi_t\right]}\\[5pt] {\mathbb{E}\big[\psi_t^2\big]}\\[5pt] {\mathbb{E}\big[\psi_t^3\big]}\end{array}\right]+\left[\begin{array}{c}\lambda {\mathbb{E}\!\left[J_1\right]} \\[5pt] \lambda {\mathbb{E}\big[J_1^2\big]} \\[5pt] \lambda {\mathbb{E}\big[J_1^3\big]}\end{array}\right].\end{equation*}

By use of the binomial theorem, we can observe that if ${\mathbb{E}\!\left[J_1^n\right]} < \infty$ then the nth moment of the shot noise process satisfies

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[\psi_t^n\right]}={\mathbb{E}\!\left[\lambda \!\left( (\psi_t + J_1)^n - \psi_t^n \right) - n \beta \psi_t^n\right]}=\sum_{k=0}^{n-1} \!\left(\begin{array}{c}n\\k\end{array}\right){\mathbb{E}\!\left[J_1^{n-k}\right]}{\mathbb{E}\!\left[\psi_t^k\right]} - n \beta {\mathbb{E}\!\left[\psi_t^n\right]},\end{equation*}

which means that the n-dimensional system is equal to the preceding one augmented below by the row vector

\begin{equation*}\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}n\lambda{\mathbb{E}\big[J_1^{n-1}\big]},& \scriptsize\left(\begin{array}{c}n\\2\end{array}\right)\!\lambda{\mathbb{E}\big[J_1^{n-2}\big]},& \scriptsize\left(\begin{array}{c}n\\3\end{array}\right)\!\lambda{\mathbb{E}\big[J_1^{n-3}\big]},& \dots,& \scriptsize\left(\begin{array}{c}n\\n - 2\end{array}\right)\!\lambda {\mathbb{E}\big[J_1^2\big]},& n\lambda{\mathbb{E}\!\left[J_1\right]} & -n\beta\end{array}\right],\end{equation*}

and to the right by zeros. Bringing this together, this now leads us to Corollary 2.

Corollary 2. Let $\psi_t$ be the intensity of a shot noise process with epochs given by a Poisson process with rate $\lambda > 0$ , jump sizes drawn from the sequence of independent and identically distributed random variables $\big\{J_i \mid i \in \mathbb{Z}^+\big\}$ , and exponential decay at rate $\beta > 0$ . If ${\mathbb{E}\!\left[J_1^n\right]} < \infty$ , the nth moment of $\psi_t$ is given by

(19) \begin{align}{\mathbb{E}\!\left[\psi_t^n\right]}&=\textbf{m}^{\psi}_{n}\!\left(\textbf{M}^{\psi}_{n-1} + n \beta\textbf{I}\right)^{-1}\!\left(e^{\textbf{M}^{\psi}_{n-1} t} - e^{-n \beta t}\textbf{I}\right)\left(\textbf{x}_{n-1}\!\left(\psi_0\right)-\frac{\lambda \textbf{j}_{n-1}}{n\beta}\right)+\psi_0^n e^{-n\beta t}\nonumber\\&\quad+\frac{\lambda{\mathbb{E}\!\left[J_1^n\right]}}{n\beta}\!\left(1-e^{-n\beta t}\right)-\frac{\lambda \textbf{m}^{\psi}_{n}}{n\beta}\!\left(\textbf{M}^{\psi}_{n-1}\right)^{-1}\Big(\textbf{I}-e^{\textbf{M}^{\psi}_{n-1}t}\Big)\textbf{j}_{n-1},\end{align}

for all $t \geq 0$ and $n \in \mathbb{Z}^+$ , where $\textbf{j}_{n} \in \mathbb{R}^n$ is such that $\left(\textbf{j}_{n}\right)_i = {\mathbb{E}\!\left[J_1^i\right]}$ ; $\textbf{M}_n^\psi \in \mathbb{R}^{n \times n}$ is recursively defined by

\begin{equation*}\textbf{M}_n^\psi=\left[\begin{array}{c@{\quad}c}\textbf{M}_{n-1}^\psi & \textbf{0}_{n-1 \times 1} \\[5pt]\textbf{m}_n^\psi & - n\beta\end{array}\right],\end{equation*}

with the row vector $\textbf{m}_n^{\psi} \in \mathbb{R}^{n-1}$ defined so that $\big(\textbf{m}_n^{\psi}\big)_i = \scriptsize\left(\begin{array}{c}n\\i\end{array}\right)\!\lambda {\mathbb{E}\big[J_1^{n-i}\big]}$ , and with $\textbf{M}_1^\psi = -\beta$ ; and $\textbf{x}_n(a) \in \mathbb{R}^n$ is such that $\left(\textbf{x}_n(a)\right)_i = a^i$ . In steady state, the $(n+1)$ th moment of the shot noise process is given by

(20) \begin{align}\lim_{t \to \infty} {\mathbb{E}\!\left[\psi_t^n\right]}&=\frac{\lambda}{n\beta}\!\left({\mathbb{E}\!\left[J_1^n\right]} - \textbf{m}_n^\psi \!\left( \textbf{M}_{n-1}^\psi\right)^{-1} \textbf{j}_{n-1} \right)\end{align}

for $n \geq 2$ , where $\lim_{t \to \infty} {\mathbb{E}\!\left[\psi_t\right]} = \frac{\lambda}{\beta}{\mathbb{E}\!\left[J_1\right]}$ . Moreover, if ${\mathbb{E}\big[J_1^{n+1}\big]} < \infty$ , the $(n+1)$ th moment of the shot noise process is given by the recursion

(21) \begin{align}\lim_{t \to \infty} {\mathbb{E}\big[\psi_t^{n+1}\big]}&=\frac{1}{(n+1)\beta}\!\left(\textbf{m}_{n+1}^\psi \textbf{s}_n^\psi + {\mathbb{E}\big[J_1^{n+1}\big]}\right),\end{align}

for all $n \in \mathbb{Z}^+$ , where $\textbf{s}_n^\psi \in \mathbb{R}^n$ is the vector of steady-state moments defined so that $\big(\textbf{s}_n^\psi\big)_i = \lim_{t \to \infty} {\mathbb{E}\!\left[\psi_t^i\right]}$ for $1 \leq i \leq n$ .

3.4. Application to Itô diffusions

For our third example, we consider an Itô diffusion; see e.g. [Reference Øksendal45] for an overview. Let $S_t$ be given by the stochastic differential equation

\begin{equation*}\textrm{d}S_t = g(S_t) \textrm{d}t + h(S_t) \textrm{d}B_t,\end{equation*}

where $B_t$ is a Brownian motion and $g({\cdot})$ and $h({\cdot})$ are real-valued functions. By Theorem 7.9 of [Reference Øksendal45], the infinitesimal generator for this process is given by

\begin{equation*}\mathcal{L}\,f(S_t)=g(S_t) f^{\prime}(S_t)+\frac{h(S_t)^2}{2} f^{\prime\prime}(S_t),\end{equation*}

where f ^′′(x) is the second derivative of $f({\cdot})$ evaluated at x. Because we will be considering functions of the form $f(x) = x^n$ for $n \in \mathbb{Z}^+$ , we will now specify the forms of $g({\cdot})$ and $h({\cdot})$ to be $g(x) = \mu + \theta x$ for some $\mu \in \mathbb{R}$ and $\theta \in \mathbb{R}$ and $h(x) = \sigma x^{\gamma/ 2}$ for some $\sigma \in \mathbb{R}$ and $\gamma \in \{0,1,2\}$ . One can note that this encompasses a myriad of relevant stochastic processes including many that are popular in the financial models literature, such as Ornstein–Uhlenbeck processes, geometric Brownian motion, and Cox–Ingersoll–Ross processes. In this case, the infinitesimal generator becomes

\begin{equation*}\mathcal{L}\,f(S_t)=\left(\mu + \theta S_t\right) f^{\prime}(S_t)+\frac{\sigma^2 S_t^{\gamma}}{2} f^{\prime\prime}(S_t),\end{equation*}

meaning that we can express the ODE for the mean as

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[S_t\right]}=\mu + \theta {\mathbb{E}\!\left[S_t\right]},\end{equation*}

and similarly the second moment will be given by the solution to

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\big[S_t^2\big]}={\mathbb{E}\big[2(\mu + \theta S_t) S_t + \sigma^2 S_t^{\gamma}\big]}=2\mu {\mathbb{E}\!\left[S_t\right]} + 2\theta {\mathbb{E}\big[S_t^2\big]} + \sigma^2 {\mathbb{E}\!\left[S_t^\gamma\right]}.\end{equation*}

For the sake of example, we now let $\gamma = 1$ , as is the case in the Cox–Ingersoll–Ross process. Then the first two transient moments of $S_t$ will be given by the solution to the system

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}\!\left[\begin{array}{c}{\mathbb{E}\!\left[S_t\right]} \\[5pt]{\mathbb{E}\big[S_t^2\big]}\\[2pt]\end{array}\right]=\left[\begin{array}{c@{\quad}c}\theta & 0 \\[5pt]2 \mu + \sigma^2 & 2\theta\end{array}\right]\left[\begin{array}{c}{\mathbb{E}\!\left[S_t\right]} \\[5pt]{\mathbb{E}\big[S_t^2\big]}\\[2pt]\end{array}\right]+\left[\begin{array}{c}\mu \\[5pt]0\end{array}\right].\end{equation*}

By observing that the third moment differential equation is

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\big[S_t^3\big]}={\mathbb{E}\Big[3(\mu + \theta S_t) S_t^2 + 3\sigma^2 S_t^{\gamma + 1}\Big]}=3\mu {\mathbb{E}\big[S_t^2\big]} + 3\theta {\mathbb{E}\big[S_t^3\big]} + 3 \sigma^2 {\mathbb{E}\big[S_t^{\gamma + 1}\big]},\end{equation*}

we can note that the third moment system for $\gamma = 1$ is

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}\left[\begin{array}{c}{\mathbb{E}\!\left[S_t\right]} \\[5pt]{\mathbb{E}\big[S_t^2\big]} \\[5pt]{\mathbb{E}\big[S_t^3\big]}\\[2pt]\end{array}\right]=\left[\begin{array}{c@{\quad}c@{\quad}c}\theta & 0 & 0 \\[5pt]2 \mu + \sigma^2 & 2\theta & 0 \\[5pt]0 & 3\mu + 3 \sigma^2 & 3\theta\end{array}\right]\left[\begin{array}{c}{\mathbb{E}\!\left[S_t\right]} \\[5pt]{\mathbb{E}\big[S_t^2\big]} \\[5pt]{\mathbb{E}\big[S_t^3\big]}\\[2pt]\end{array}\right]+\left[\begin{array}{c}\mu \\[5pt]0 \\[5pt]0\end{array}\right],\end{equation*}

and this showcases the matryoshkan nesting structure, as the second moment system is contained within the third. Because the general nth moment for $n \geq 2$ has differential equation given by

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[S_t^n\right]}&={\mathbb{E}\!\left[n(\mu + \theta S_t) S_t^{n-1} + \frac{n(n-1)\sigma^2}{2} S_t^{n + \gamma - 2}\right]}\\&=n\mu {\mathbb{E}\big[S_t^{n-1}\big]} + n\theta {\mathbb{E}\!\left[S_t^n\right]} + \frac{n(n-1)\sigma^2}{2} {\mathbb{E}\Big[S_t^{n + \gamma - 2}\Big]},\end{align*}

we can see that the $(n-1)$ th system can be augmented below by the row vector for $\gamma = 1$ ,

\begin{equation*}\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}0, & 0, & \dots, & 0, & n\mu + \frac{n(n-1)\sigma^2}{2}, & n\theta\end{array}\right],\end{equation*}

and to the right by zeros. Through this observation, we can now give the moments of Itô diffusions in Corollary 3.

Corollary 3. Let $S_t$ be an Itô diffusion that satisfies the stochastic differential equation

(22) \begin{align}\textrm{d}S_t=(\mu + \theta S_t)\textrm{d}t+\sigma S_t^{\gamma / 2} \textrm{d} B_t,\end{align}

where $B_t$ is a Brownian motion, $\mu, \theta, \sigma \in \mathbb{R}$ , and $\gamma \in \{0, 1, 2\}$ . Then the nth moment of $S_t$ is given by

(23) \begin{align}{\mathbb{E}\!\left[S_t^n\right]}&=\textbf{m}^S_{n}\Big(\textbf{M}^S_{n-1} - \chi_n\textbf{I}\Big)^{-1}\Big(e^{\textbf{M}^S_{n-1} t} - e^{\chi_n t}\textbf{I}\Big)\left(\textbf{x}_{n-1}(S_0)+\frac{\mu \textbf{v}_1 + \sigma^2\mathbb{I}_{\{\gamma = 0\}}\textbf{v}_2}{\chi_n}\right)\nonumber\\&\quad+S_0^n e^{\chi_n t}-\big(\mu \mathbb{I}_{\{n=1\}} + \sigma^2\mathbb{I}_{\{\gamma = 0, n= 2\}}\big)\frac{1-e^{\chi_n t}}{\chi_n}\nonumber\\&\quad+\frac{\textbf{m}^S_{n}}{\chi_{n}}\!\left(\textbf{M}^S_{n-1} \right)^{-1}\!\left(\textbf{I}-e^{\textbf{M}^S_{n-1}t}\right)\left(\mu\textbf{v}_1 + \sigma^2\mathbb{I}_{\{\gamma = 0\}}\textbf{v}_2\right),\end{align}

for all $t \geq 0$ and $n \in \mathbb{Z}^+$ , where $\chi_n = n\theta + \frac{n}{2}(n-1)\sigma^2\mathbb{I}_{\{\gamma = 2\}}$ ;

\begin{equation*}\textbf{M}^S_n = \theta \textbf{diag}(1\,:\,n) + \mu \textbf{diag}(2\,:\,n,-1) + \frac{\sigma^2}{2}\textbf{diag}\!\left(\textbf{d}^{2-\gamma}_{n+\gamma-2}, \gamma - 2 \right)\end{equation*}

for $\textbf{d}^j_k \in \mathbb{R}^k$ such that $\big(\textbf{d}_k^j\big)_i = (\,j+i)(\,j+i-1)$ ;

\begin{equation*}\textbf{m}_n^{S}=\left[\big(\textbf{M}_n^S\big)_{n,1},\,\dots,\,\big(\textbf{M}_n^S\big)_{n,n-1}\right]\end{equation*}

is such that

\begin{equation*}\left(\textbf{m}_n^{S}\right)_{j}=\begin{cases}n\mu + \frac{n(n-1)\sigma^2}{2}\mathbb{I}_{\{\gamma = 1\}}, & j = n -1, \\[5pt]\frac{n(n-1)\sigma^2}{2}\mathbb{I}_{\{\gamma = 0\}}, & j = n - 2, \\[4pt]0, & 1 \leq j < n - 2;\end{cases}\end{equation*}

and $\textbf{x}_n(a) \in \mathbb{R}^n$ is such that $\left(\textbf{x}_n(a)\right)_i = a^i$ . If $\theta < 0$ and $\gamma \in \{0,1\}$ , then the nth steady-state moment of $S_t$ is given by

(24) \begin{align}\lim_{t \to \infty}{\mathbb{E}\!\left[S_t^n\right]}=\frac{\mu}{\chi_n}\textbf{m}_n^S \!\left(\textbf{M}_{n-1}^S\right)^{-1} \textbf{v}_1,\end{align}

for $n \geq 2$ , with $\lim_{t \to \infty} {\mathbb{E}\!\left[S_t\right]} = -\frac{\mu}{\theta}$ . Moreover, the $(n+1)$ th steady-state moment of $S_t$ is given by the recursion

(25) \begin{align}\lim_{t \to \infty} {\mathbb{E}\big[S_t^{n+1}\big]}&=-\frac{1}{\chi_n}\textbf{m}_{n+1}^S \textbf{s}_n^S,\end{align}

for all $n \in \mathbb{Z}^+$ , where $\textbf{s}_n^S \in \mathbb{R}^n$ is the vector of steady-state moments defined so that $\left(\textbf{s}_n^S\right)_i = \lim_{t \to \infty}{\mathbb{E}\!\left[S_t^i\right]}$ for $1 \leq i \leq n$ .

As a consequence of these expressions we can also gain insight on the moments of an Itô diffusion in the case of non-integer $\gamma \in [0,2]$ , as is used in volatility models such as the CEV model and the SABR model; see e.g. [Reference Henry-Labordère34]. This can be achieved through bounding the differential equations, as the nth moment of such a diffusion is again given by

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[S_t^n\right]}=n\mu {\mathbb{E}\big[S_t^{n-1}\big]} + n\theta {\mathbb{E}\!\left[S_t^n\right]} + \frac{n(n-1)\sigma^2}{2} {\mathbb{E}\Big[S_t^{n + \gamma - 2}\Big]},\end{equation*}

the rightmost term of which can be bounded above and below by

\begin{equation*} {\mathbb{E}\Big[S_t^{n + \lfloor\gamma\rfloor - 2}\Big]}\leq {\mathbb{E}\Big[S_t^{n + \gamma - 2}\Big]}\leq {\mathbb{E}\Big[S_t^{n + \lceil\gamma\rceil - 2}\Big]},\end{equation*}

and the differential equations given by substituting these bounded terms form a closed system solvable by Corollary 3. Assuming the true differential equation and the upper and lower bounds all share an initial value, the solution to the bounded equations bounds the solution to the true moment equation; see [Reference Hale and Lunel32].

3.5. Application to growth–collapse processes

For a fourth example, we consider growth–collapse processes with Poisson-driven shocks. These processes have been studied in a variety of contexts; see e.g. [Reference Boxma, Perry, Stadje and Zacks7, Reference Kella38, Reference Kella and Löpker39, Reference Boxma, Kella and Perry6]. More recently, these processes and their related extensions have seen renewed interest in connection with the study of the crypto-currency Bitcoin; see for example [Reference Frolkova and Mandjes28, Reference Koops40, Reference Javier and Fralix37, Reference Fralix27]. While growth–collapse processes can be defined in many different ways, for this example we use a definition in the style of Section 4 from [Reference Boxma, Perry, Stadje and Zacks7]. We let $Y_t$ be the state of the growth–collapse model and let $\{U_i \mid i \in \mathbb{Z}^+\}$ be a sequence of independent $\textrm{Uni}(0,1)$ random variables that are also independent from the state and history of the growth–collapse process. From Proposition 1 of [Reference Boxma, Perry, Stadje and Zacks7], the infinitesimal generator of $Y_t$ is given by

\begin{equation*}\mathcal{L}\,f(Y_t)= \lambda f^{\prime}(Y_t) + \mu \int_0^1 \!\left( f( u Y_t ) - f( Y_t ) \right) \textrm{d}u.\end{equation*}

Thus, $Y_t$ experiences linear growth at rate $\lambda > 0$ throughout time, but it also collapses at epochs given by a Poisson process with rate $\mu > 0$ . At the ith collapse epoch the process falls to a fraction of its current level; specifically it jumps down to $U_i Y_t$ . Using the infinitesimal generator, we can see that the mean of this growth–collapse process satisfies

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[Y_t\right]}=\lambda + \mu \!\left( {\mathbb{E}\!\left[U_1 Y_t\right]} - {\mathbb{E}\!\left[Y_t\right]} \right)=\lambda - \frac{\mu}{2} {\mathbb{E}\!\left[Y_t\right]},\end{equation*}

and its second moment satisfies

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\big[Y_t^2\big]}=2 \lambda {\mathbb{E}\!\left[Y_t\right]} + \mu \!\left( {\mathbb{E}\big[U_1^2 Y_t^2\big]} - {\mathbb{E}\big[Y_t^2\big]} \right)=2\lambda {\mathbb{E}\!\left[Y_t\right]} - \frac{2 \mu}{3} {\mathbb{E}\big[Y_t^2\big]}.\end{equation*}

Therefore, we can write the linear system of differential equations for the second moment of this growth–collapse process as

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}\!\left[\begin{array}{c}{\mathbb{E}\!\left[Y_t\right]}\\[5pt] {\mathbb{E}\big[Y_t^2\big]}\\[2pt]\end{array}\right]=\left[\begin{array}{c@{\quad}c}-\frac{\mu}2 & 0 \\[5pt] 2\lambda & - \frac{2\mu}{3}\\[2pt]\end{array}\right]\left[\begin{array}{c}{\mathbb{E}\!\left[Y_t\right]}\\[5pt] {\mathbb{E}\big[Y_t^2\big]}\\[2pt]\end{array}\right]+\left[\begin{array}{c}\lambda \\[5pt] 0\end{array}\right].\end{equation*}

Moving to the third moment, via the infinitesimal generator we write its differential equation as

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\big[Y_t^3\big]}=3 \lambda {\mathbb{E}\big[Y_t^2\big]} + \mu \!\left( {\mathbb{E}\big[U_1^3 Y_t^3\big]} - {\mathbb{E}\big[Y_t^3\big]} \right)=3\lambda {\mathbb{E}\big[Y_t^2\big]} - \frac{3 \mu}{4} {\mathbb{E}\big[Y_t^3\big]},\end{equation*}

which shows that the system of differential equations for the third moment is

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}\!\left[\begin{array}{c}{\mathbb{E}\!\left[Y_t\right]}\\[5pt] {\mathbb{E}\big[Y_t^2\big]}\\[5pt] {\mathbb{E}\big[Y_t^3\big]}\\[2pt]\end{array}\right]= \left[\begin{array}{c@{\quad}c@{\quad}c}-\frac{\mu}2 & 0 & 0\\[5pt] 2\lambda & - \frac{2\mu}{3} & 0\\[5pt] 0 & 3\lambda & - \frac{3\mu}{4}\\[2pt]\end{array}\right] \left[\begin{array}{c}{\mathbb{E}\!\left[Y_t\right]}\\[5pt] {\mathbb{E}\big[Y_t^2\big]}\\[5pt] {\mathbb{E}\big[Y_t^3\big]}\\[2pt]\end{array}\right]+ \left[\begin{array}{c}\lambda \\[5pt] 0\\[5pt] 0\end{array}\right];\end{equation*}

this obviously encapsulates the system for the first two moments. We can note that the general nth moment will satisfy

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[Y_t^n\right]}=n \lambda {\mathbb{E}\big[Y_t^{n-1}\big]} + \mu \!\left( {\mathbb{E}\!\left[U_1^n Y_t^n\right]} - {\mathbb{E}\!\left[Y_t^n\right]} \right)=n\lambda {\mathbb{E}\big[Y_t^{n-1}\big]} - \frac{n \mu}{n+1} {\mathbb{E}\!\left[Y_t^n\right]},\end{equation*}

and thus the system for the nth moment is given by appending the row vector

\begin{equation*}\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}0, &, 0, & \dots, & 0, & n\lambda, & -\frac{n\mu}{n+1}\end{array}\right]\end{equation*}

below the matrix from the $(n-1)$ th system augmented by zeros on the right. Following this derivation, we reach the general expressions for the moments in Corollary 4. Furthermore, we can note that because of the relative simplicity of this particular structure, we are able to solve the recursion for the steady-state moments and give these terms explicitly.

Corollary 4. Let $Y_t$ be a growth–collapse process with growth rate $\lambda > 0$ and uniformly sized collapses occurring according to a Poisson process with rate $\mu > 0$ . Then the nth moment of $Y_t$ is given by

(26) \begin{align}{\mathbb{E}\!\left[Y_t^n\right]}&=n\lambda\textbf{v}_{n-1}^{\textrm{T}}\!\left(\textbf{M}^Y_{n-1} + \frac{n\mu}{n+1}\textbf{I}\right)^{-1}\!\left(e^{\textbf{M}^Y_{n-1} t} - e^{-\frac{n\mu t}{n+1}}\textbf{I}\right)\left(\textbf{x}_{n-1}(y_0)-\frac{(n+1)\lambda \textbf{v}_1}{n\mu}\right)\nonumber\\&\quad+y_0^n e^{-\frac{n\mu t}{n+1}}+\frac{(n+1)\lambda \mathbb{I}_{\{n=1\}}}{n\mu}\!\left(1-e^{-\frac{n\mu t}{n+1}}\right)\nonumber\\&\quad-\frac{(n+1)\lambda^2}{\mu}\textbf{v}^{\textrm{T}}_{n-1}\!\left(\textbf{M}^Y_{n-1}\right)^{-1}\!\left(\textbf{I}-e^{\textbf{M}^Y_{n-1}t}\right)\textbf{v}_1,\end{align}

where $y_0$ is the initial value of $Y_t$ , $\textbf{x}_n(a) \in \mathbb{R}^n$ is such that $\left(\textbf{x}_n(a)\right)_i = a^i$ ,

\begin{equation*}\textbf{M}_n^Y = \lambda \textbf{diag}\!\left(2\,:\,n,-1\right) - \mu \textbf{diag}\!\left(\frac{1}{2} \,:\, \frac{n}{n+1}\right),\end{equation*}

and

\begin{equation*}\textbf{m}_n^Y = \left[\left(\textbf{M}_n^Y\right)_{n,1}, \dots, \left(\textbf{M}_n^Y\right)_{n,n-1} \right]\end{equation*}

is such that $\textbf{m}_n^Y = n\lambda\textbf{v}_{n-1}^{\textrm{T}}$ . Moreover, the nth steady-state moment of $Y_t$ is given by

(27) \begin{align}\lim_{t \to \infty}{\mathbb{E}\!\left[Y_t^n\right]}&=(n+1)!\left(\frac{\lambda}{\mu}\right)^n,\end{align}

for $n \in \mathbb{Z}^+$ .

3.6. Application to ephemerally self-exciting processes

As a final detailed example of the applicability of matryoshkan matrices, we now consider a stochastic process we have analyzed in [Reference Daw and Pender19]. This process is a linear birth–death–immigration process, which has been of interest in classical teletraffic theory; see for example [Reference Wallström46], [Reference Iversen35], and Section 8.4 of [Reference Iversen36]. By comparison, our interests in [Reference Daw and Pender19] are in alternate representations of self-exciting processes and alternate constructions of the original self-exciting process, the Hawkes process. Here, the occurrence of an arrival increases the arrival rate by an amount $\alpha > 0$ , as in the Hawkes process, and this increase expires after an exponentially distributed duration with some rate $\mu > \alpha$ . Accordingly, this process is an ephemerally self-exciting process. Given a baseline intensity $\nu^* > 0$ , let $Q_t$ be such that new arrivals occur at rate $\nu^* + \alpha Q_t$ and then the overall rate until the next excitement expiration is $\mu Q_t$ . One can then think of $Q_t$ as the number of entities still causing active excitement at time $t \geq 0$ . We will refer to $Q_t$ as the number in system for this ephemerally self-exciting process. The infinitesimal generator for a function $f \,:\, \mathbb{N} \to \mathbb{R}$ is thus

\begin{equation*}\mathcal{L}\,f(Q_t) = \left(\nu^* + \alpha Q_t\right) \left(\,f(Q_t + 1) - f(Q_t)\right) + \mu Q_t \!\left(\,f(Q_t - 1) - f(Q_t)\right),\end{equation*}

which again captures the dynamics of the process, as the first term on the right-hand side is the product of the up-jump rate and the change in function value upon an increase in the process, while the second term is the product of the down-jump rate and the corresponding process decrease. This now yields an ODE for the mean given by

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[Q_t\right]}=\nu^* + \alpha{\mathbb{E}\!\left[Q_t\right]} - \mu {\mathbb{E}\!\left[Q_t\right]}=\nu^*-(\mu - \alpha) {\mathbb{E}\!\left[Q_t\right]},\end{equation*}

while the second moment will satisfy

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\big[Q_t^2\big]}&={\mathbb{E}\!\left[\left(\nu^* + \alpha Q_t\right)\left((Q_t+1)^2 - Q_t^2\right) + \mu Q_t \!\left((Q_t-1)^2 - Q_t^2\right)\right]}\\&=\left(2\nu^* + \mu + \alpha\right){\mathbb{E}\!\left[Q_t\right]}+ \nu^*- 2(\mu - \alpha){\mathbb{E}\big[Q_t^2\big]}.\end{align*}

Thus, the first two moments are given by the solution to the linear system

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}\!\left[\begin{array}{c}{\mathbb{E}\!\left[Q_t\right]}\\[4pt] {\mathbb{E}\big[Q_t^2\big]}\\[2pt]\end{array}\right]=\left[\begin{array}{c@{\quad}c}-(\mu - \alpha) & 0 \\[4pt] 2\nu^* + \mu + \alpha & -2(\mu-\alpha)\end{array}\right]\left[\begin{array}{c}{\mathbb{E}\!\left[Q_t\right]}\\[4pt] {\mathbb{E}\big[Q_t^2\big]}\\[2pt]\end{array}\right]+\left[\begin{array}{c}\nu^* \\[4pt] \nu^*\end{array}\right],\end{equation*}

and by observing that the third moment differential equation is

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\big[Q_t^3\big]}&={\mathbb{E}\!\left[\left(\nu^* + \alpha Q_t\right)\big((Q_t+1)^3 - Q_t^3\big) + \mu Q_t \big((Q_t-1)^3 - Q_t^3\big)\right]}\\&=\left(3\nu^* + 3\alpha + 3\mu\right){\mathbb{E}\big[Q_t^2\big]}+\left(3\nu^* + \alpha - \mu\right){\mathbb{E}\!\left[Q_t\right]}+\nu^*- 3(\mu - \alpha){\mathbb{E}\big[Q_t^3\big]},\end{align*}

we can see that the third moment system does indeed encapsulate the second moment system:

\begin{equation*}\frac{\textrm{d}}{\textrm{d}t}\!\left[\begin{array}{c}{\mathbb{E}\!\left[Q_t\right]}\\[5pt] {\mathbb{E}\big[Q_t^2\big]} \\[5pt] {\mathbb{E}\big[Q_t^3\big]}\\[2pt]\end{array}\right]=\left[\begin{array}{c@{\quad}c@{\quad}c}-(\mu - \alpha) & 0 & 0 \\[5pt] 2\nu^* + \mu + \alpha & -2(\mu-\alpha) & 0 \\[5pt] 3\nu^* + \alpha - \mu & 3\nu^* + 3\alpha + 3\mu & - 3(\mu-\alpha)\\[2pt]\end{array}\right] \left[\begin{array}{c} {\mathbb{E}\!\left[Q_t\right]}\\[5pt] {\mathbb{E}\big[Q_t^2\big]} \\[5pt] {\mathbb{E}\big[Q_t^3\big]}\\[2pt]\end{array}\right]+ \left[\begin{array}{c}\nu^* \\[5pt] \nu^* \\[5pt] \nu^*\end{array}\right].\end{equation*}

In general, the nth moment is given by the solution to

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[Q_t^n\right]}&={\mathbb{E}\!\left[(\nu^* + \alpha Q_t) \left((Q_t + 1)^n - Q_t^n\right) + \mu Q_t \!\left((Q_t - 1)^n - Q_t^n\right)\right]}\\&=\nu^*+\nu^* \sum_{k=1}^{n-1}\!\left(\begin{array}{c}n\\k\end{array}\right) {\mathbb{E}\!\left[Q_t^k\right]}+\alpha \sum_{k=1}^{n}\!\left(\begin{array}{c}n\\k - 1\end{array}\right) {\mathbb{E}\!\left[Q_t^k\right]}\\& \quad+\mu \sum_{k=1}^{n} \!\left(\begin{array}{c}n\\k - 1\end{array}\right) {\mathbb{E}\!\left[Q_t^{k}\right]}({-}1)^{n-k-1},\end{align*}

which means that the nth system is given by augmenting the previous system below by

\begin{equation*}\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}n\nu^* + \alpha + \mu({-}1)^n,&\!\!\scriptsize\left(\begin{array}{c}n\\2\end{array}\right)\!\nu^* + n\alpha + n\mu({-}1)^{n-1},&\!\!\dots,&\!\! n\nu^* + \scriptsize\left(\begin{array}{c}n\\2\end{array}\right)\!\alpha + \scriptsize\left(\begin{array}{c}n\\2\end{array}\right)\!\mu,&\!\! -n(\mu-\alpha)\end{array}\right],\end{equation*}

and to the right by zeros. By comparing this row vector to the definition of the matryoshkan Pascal matrices in Equation (15), we arrive at explicit forms for the moments of this process, shown now in Corollary 5.

Corollary 5. Let $Q_t$ be the number in system for an ephemerally self-exciting process with baseline intensity $\nu^* > 0$ , intensity jump $\alpha > 0$ , and duration rate $\mu > \alpha$ . Then the nth moment of $Q_t$ is given by

(28) \begin{align}{\mathbb{E}\!\left[Q_t^n\right]}&=\textbf{m}^Q_{n}\!\left(\textbf{M}^Q_{n-1} + n(\mu-\alpha)\textbf{I}\right)^{-1}\!\left(e^{\textbf{M}^Q_{n-1} t} - e^{-n(\mu-\alpha)t}\textbf{I}\right)\left(\textbf{x}_{n-1}(Q_0)-\frac{\nu^* \textbf{v}}{n(\mu-\alpha)}\right)\nonumber\\&\quad+Q_0^n e^{-n(\mu-\alpha)t}+\frac{\nu^*}{n(\mu-\alpha)}\!\left(1-e^{-n(\mu-\alpha)t}\right)\nonumber\\&\quad-\frac{\nu^* \textbf{m}^Q_{n}}{n(\mu-\alpha)}\!\left(\textbf{M}^Q_{n-1}\right)^{-1}\!\left(\textbf{I}-e^{\textbf{M}^Q_{n-1}t}\right)\textbf{v},\end{align}

for all $t \geq 0$ and $n \in \mathbb{Z}^+$ , where $\textbf{M}^Q_n = \nu^*\boldsymbol{\mathcal{P}}_n(1)\textbf{diag}(\textbf{v},-1) + \alpha\boldsymbol{\mathcal{P}}_n(1) + \mu\boldsymbol{\mathcal{P}}_n({-}1)$ ,

\begin{equation*}\textbf{m}_n^{Q}=\left[\left(\textbf{M}_n^Q\right)_{n,1},\,\dots,\,\left(\textbf{M}_n^Q\right)_{n,n-1}\right]\end{equation*}

is such that

\begin{equation*}\left(\textbf{m}_n^{Q}\right)_{j}=\left(\begin{array}{c}n\\j\end{array}\right) \nu^*+\left(\begin{array}{c}n\\j - 1\end{array}\right)\alpha+\left(\begin{array}{c}n\\j - 1\end{array}\right)\mu ({-}1)^{n-j-1},\end{equation*}

and $\textbf{x}_n(a) \in \mathbb{R}^n$ is such that $\left(\textbf{x}_n(a)\right)_i = a^i$ . In steady state, the nth moment of $Q_t$ is given by

(29) \begin{align}\lim_{t \to \infty}{\mathbb{E}\!\left[Q_t^n\right]}=\frac{\nu^*}{\mu - \alpha}\!\left(1 - \textbf{m}_n^Q \!\left(\textbf{M}_{n-1}^Q\right)^{-1} \textbf{v}\right),\end{align}

for $n \geq 2$ , with $\lim_{t \to \infty} {\mathbb{E}\!\left[Q_t\right]} = \frac{\nu^*}{\mu-\alpha}$ . Moreover, the $(n+1)$ th steady-state moment of the ephemerally self-exciting process is given by the recursion

(30) \begin{align}\lim_{t \to \infty} {\mathbb{E}\big[Q_t^{n+1}\big]}&=\frac{1}{(n+1)(\mu-\alpha)}\Big(\textbf{m}_{n+1}^Q \textbf{s}_n^Q + \nu^*\Big),\end{align}

for all $n \in \mathbb{Z}^+$ , where $\textbf{s}_n^Q \in \mathbb{R}^n$ is the vector of steady-state moments defined so that $\big(\textbf{s}_n^Q\big)_i = \lim_{t \to \infty}{\mathbb{E}\!\left[Q_t^i\right]}$ for $1 \leq i \leq n$ .

Let us note that the ephemerally self-exciting process is known to be negative binomially distributed in steady state [Reference Daw and Pender19]. Hence, while the vector and recursive computations in Equations (29) and (30), respectively, may be useful for higher moments, for lower moments there may be closed-form expressions available.

3.7. Additional applications by combination and permutation

While the preceding examples are the the only detailed examples we include in this paper, we can note that these matryoshkan matrix methods can be applied to many other settings. In fact, one can observe that these example derivations can be applied directly to processes that feature a combination of their structures, such as the dynamic contagion process introduced in [Reference Dassios and Zhao17]. The dynamic contagion process is a point process that is both self-excited and externally excited, meaning that its intensity experiences jumps driven both by its own activity and by the activity of an exogenous Poisson process. In this way, the process combines the behavior of the Hawkes and shot noise processes. Hence, its infinitesimal generator can be written using a combination of expressions used in Subsections 3.2 and 3.3, implying that all moments of the process can be calculated through this methodology. Similarly, these methods can also be readily applied to processes that combine dynamics from Hawkes processes and from Itô diffusions, such as affine point processes. These processes, studied in e.g. [Reference Errais, Giesecke and Goldberg22, Reference Zhang, Blanchet, Giesecke and Glynn48, Reference Gao and Zhu30], feature both self-excitement and diffusion behavior and thus have an infinitesimal generator that can be expressed using terms from the generators for Hawkes and Itô processes. One could also study the combination of externally driven jumps and diffusive behavior, such as in affine jump diffusions; see e.g. [Reference Duffie, Pan and Singleton20]. Of course, one can also consider permutations of the model features seen in our examples, such as trading fixed-size jumps for random ones to form marked Hawkes processes or changing to randomly sized batches of arrivals in the ephemerally self-exciting process. In general, the key requirement from the assumptions in Theorem 1 is the closure of the system of moment differential equations specified in Equation (7). This is equivalent to having the infinitesimal generator of any polynomial being a polynomial of order no more than the original, which again connects this work to the literature on polynomial processes. In fact, the notion of closure of polynomial processes under combinations of infinitesimal generators has been studied by [Reference Cuchiero, Larsson and Svaluto-Ferro12]. In summary, infinitesimal generators of the form

\begin{align*}\mathcal{L}\,f(X_t)&=\underbrace{(\alpha_0 + \alpha_1 X_t ) \int_0^\infty \!\left(\,f(X_t + a) - f(X_t) \right)\! \textrm{d}F_{A}(a) }_{\text{Up-jumps}}+\underbrace{(\alpha_2 + \alpha_3 X_t ) f^{\prime}(X_t)}_{\text{Drift, decay, or growth}}\\&\quad+\underbrace{(\alpha_4 + \alpha_5 X_t ) \int_0^\infty \!\left(\,f(X_t - b) - f(X_t) \right)\! \textrm{d}F_{B}(b)}_{\text{Down-jumps}}+\underbrace{\big(\alpha_6 + \alpha_7 X_t + \alpha_8 X_t^2 \big) f^{\prime\prime}(X_t)}_{\text{Diffusion}}\\&\quad+\underbrace{\alpha_9 \int_0^\infty \!\left( f(c X_t) - f(X_t) \right)\! \textrm{d}F_C(c)}_{\text{Expansion or collapse}}\end{align*}

can be handled by this methodology, where $\alpha_j \in \mathbb{R}$ for all j and where the sequences $\{A_i\}$ , $\{B_i\}$ , and $\{C_i\}$ are of mutually independent random variables with respective cumulative distribution functions $F_A({\cdot})$ , $F_B({\cdot})$ , and $F_C({\cdot})$ . Finally we note that this example generator need not be exclusive, as it is possible that other dynamics may also meet the closure requirements in Equation (7).