(I cross-posted this question on MO: https://mathoverflow.net/questions/450064/characteristic-function-of-stochastic-integral-of-a-pure-jump-l%c3%a9vy-process-with)
Let $X(t)$ and $Y(t)$ be two pure jump Lévy processes, so we know their characteristic functions through the Lévy-Kintchine formula as $\phi_{X(t)}(u) = \mathrm{e}^{t\psi_{X}(u)}$ and $\phi_{Y(t)}(u) = \mathrm{e}^{t\psi_{Y}(u)}$ where $\psi_{X}$ and $\psi_{Y}$ are the characteristic exponents of $X(1)$ and $Y(1)$ respectively.
I want to calculate (or approximate) the probability density function of the following stochastic integral
$$Z(t) = Z_0 + \int_{0}^t Y(s_{-}) dX(s)$$
Which I think I can obtain through calculation (and subsequent numerical inversion) of its characteristic function
$$\phi_{Z(t)(u)} = \mathbb{E}\left[\exp \left(iu \int_{0}^t Y(s_{-}) dX(s) \right) \right]$$
Characteristic function of stochastic integral of Lévy integrand and Lévy integrator
I know that for deterministic functions $f(t)$ we have (see Tankov 2003 Lemma 15.1) $$\mathbb{E}\left[\exp \left(iu \int_{0}^t f(s) dX(s) \right) \right] = \exp \left(\int_{0}^t \psi_{X}(uf(s))ds \right)$$ So I was looking to see if it was possible to find something analogous when the integrand is a stochastic process.
Starting from writing the integral from its definition, assuming we have our grid in $t_k$, we can write the characteristic function as (I'm omitting notation of the limit for when the grid becomes continuous): \begin{align*} \phi_{Z(t)(u)} &= \mathbb{E}\left[\exp \left(iu \sum_k Y(t_k) \left[ X(t_{k+1}) - X(t_k)\right] \right) \right] \\&= \mathbb{E}\left[\exp \left(iu \sum_k \left[Y(t_k) - Y_0 + Y_0 \right] X(t_{k+1} - t_k) \right) \right] \\&= \mathbb{E}\left[\exp \left(iu \sum_k \left[Y(t_k) - Y_0 \right] X(t_{k+1} - t_k) + iu\sum_k Y_0 X(t_{k+1} - t_k) \right) \right] \\&= \mathbb{E}\left[\exp \left(iu \sum_k \left[Y(t_k) - Y_0 \right] X(t_{k+1} - t_k)\right) \exp\left(iu\sum_k Y_0 X(t_{k+1} - t_k) \right) \right] \\&= \mathbb{E}\left[ \prod_k \exp \left( iu \left[Y(t_k) - Y_0\right] X(t_{k+1} - t_k) \right) \prod_k \exp \left(iu Y_0 X(t_{k+1} - t_k) \right) \right] \end{align*}
This looks like the characteristic function of the sum of two stochastic integrals, which in hindsight seems obvious as this manipulation could have been done in the stochastic integral itself before writing the expression for its characteristic function, and I think I can separate the expectations so it continues as
\begin{align} \begin{split} & \phi_{Z(t)(u)} = \mathbb{E}\left[ \prod_k \exp \left( iu \left[Y(t_k) - Y_0\right] X(t_{k+1} - t_k) \right)\right] \mathbb{E}\left[\prod_k \exp \left(iu Y_0 X(t_{k+1} - t_k) \right) \right] \end{split} \end{align}
I then take the products outside the expectations because of the independence of increments of Lévy processes, and proceed to use the formula for the characteristic function of a product of random variables
\begin{align} \phi_{Z(t)(u)} & = \prod_k \mathbb{E}\left[ \exp \left( iu \left[Y(t_k) - Y_0\right] X(t_{k+1} - t_k) \right) \right]\prod_k \mathbb{E}\left[ \exp \left( iu Y_0 X(t_{k+1} - t_k) \right) \right] \\ & = \prod_k \mathbb{E}\left[\exp \left( \psi_{X(t_{k+1}-t_k)}\left(u \left[Y(t_k) - Y_0\right] \right) \right) \right] \prod_k \mathbb{E}\left[\exp \left( \psi_{X(t_{k+1}-t_k)}\left(u Y_0\right) \right) \right] \\ & = \prod_k \mathbb{E}\left[\exp \left( (t_{k+1}-t_k)\psi_{X}\left(u \left[Y(t_k) - Y_0\right] \right) \right) \right] \prod_k \mathbb{E}\left[\exp \left( (t_{k+1}-t_k)\psi_{X}\left(u Y_0\right) \right) \right] \\ & = \mathbb{E}\left[\exp \left( \sum_k (t_{k+1}-t_k)\psi_{X}\left(u \left[Y(t_k) - Y_0\right] \right) \right) \right] \mathbb{E}\left[\exp \left( \sum_k (t_{k+1}-t_k)\psi_{X}\left(u Y_0 \right) \right) \right] \\ & \rightarrow \mathbb{E}\left[\exp \left( \int_0^t \psi_{X}\left(u \left[Y(s) - Y_0\right]\right) ds \right) \right] \mathbb{E}\left[\exp \left( \int_0^t \psi_{X}\left(u Y_0\right) ds \right) \right] \\ & = \mathbb{E}\left[\exp \left( \int_0^t \psi_{X}\left(u \left[Y(s) - Y_0\right]\right) ds \right) \right] \exp \left( \int_0^t \psi_{X}\left(u Y_0\right) ds \right) \end{align}
Now we need to work on the expectation
$$\mathbb{E}\left[\exp \left( \int_0^t \psi_{X}\left(u Y(s) - uY_0\right) ds \right) \right]$$
I will simplify by setting $Y_0 = 0$, but I believe all the following calculations will be analogous, and we would just need to replace "$ux$" by "$ux - uY_0$", "$\psi_X(0)$" by "$\psi_X(-uY_0)$" and "$Y(s)$" by "$Y(s) - Y_0$" whenever appropriate in order to recover the more general case.
Itô's formula on the integrand inside the exponent
Now I'm trying to apply Itô's formula for pure jump processes (Applebaum 2009, Lemma 4.4.5) on the integrand inside the previous expression, namely $f(t,Y(t)) = \psi_{X}( u Y(t))$. This yields something like
\begin{aligned} \begin{split} \psi_{X}( u Y(t)) & = \psi_{X}(0) + \int_0^t \frac{\partial \psi_{X}(ux)}{\partial s} \vert_{x=Y(s)} ds \\ & + \int_0^t \int_{\mathbb{R}} \left \{ \psi_{X}\left(uY(s_{-}) + uy\right) - \psi_{X}\left(u Y(s_{-})\right)\right\} N_Y(dy,ds) \end{split} \end{aligned}
I think the integral with respect to time disappears because $\psi_{X}(u) = \psi_{X(1)}(u)$ has no explicit dependency on $t$.
If I refer to the new stochastic process being studied as $A(t) := \psi_{X}( u Y(t))$, we can now write that what we want to study is its time integral $\int_0^t A(s) ds$. So there's a given area under each trajectory of the stochastic process $A(t)$ across the time interval we are integrating on, but I'm not sure how to calculate this area for a given trajectory. I was looking at a version of Fubini's theorem for stochastic integrals (Protter 2004 Theorem 64) and I found something that led me to hypothesize that I might be able to apply it in this way
$$\int_0^t A(s)ds = \int_0^t \left(\int_0^s dA(u)\right) ds = \int_0^t \left( \int_0^s du \right) dA(s) = \int_0^t s \cdot dA(s) $$
But I have no idea if this makes any sense. If it does, I can now take what we've previously calculated with the Itô formula for $dA(t)$ and write
\begin{aligned} \begin{split} & \int_0^t \psi_{X}( u Y(s))ds = \int_0^t s \cdot d\left(\psi_{X}( u Y(s))\right) \\ & = \int_0^t \int_{\mathbb{R}} s \left [ \psi_{X}\left(uY(s_{-}) + uy\right) - \psi_{X}\left(u Y(s_{-})\right)\right ] N_Y(dy,ds) \end{split} \end{aligned}
And I'm now at the step of figuring out a way to find the characteristic function of this integral w.r.t. to the Poisson random measure of $Y$. I have created a new question that pertains to the particular difficulty of not being able to eliminate the dependency on $Y(s_{-})$: Characteristic function of Itô's formula for jump processes.