2022 Prospects for low-frequency radio astronomy in S. A.

# Search for extraplanar radio emission in a sample of edge-on starburst
galaxies

C. A. Galante  11affiliation: Facultad de Ciencias Astronómicas y Geofísicas,
UNLP, Paseo del Bosque S/N, 1900, La Plata, Buenos Aires, Argentina.
22affiliation: Instituto Argentino de Radioastronomía, Berazategui, Buenos
Aires, Argentina. J. Saponara  22affiliationmark:  G. E. Romero
11affiliationmark:  22affiliationmark:  and P. Benaglia22affiliationmark:

###### Abstract

The intense star formation in starburst galaxies results in the formation of
intense winds that sweep matter up to several kpc out of the galactic plane.
These winds can be detected in many bands of the electromagnetic spectrum. In
this opportunity, we present our project to study three edge-on galaxies,
candidates to show extraplanar emission at radio wavelengths. We will describe
the criteria followed in selecting the sample and the software developed to
detect the presence of the wind. Besides, we will discuss the next steps of
this project.

La intensa formación estelar que presentan las galaxias starburst resulta en
la formación de fuertes vientos que arrastran materia hasta varios kpc por
fuera del plano galáctico, emitiendo radiación en varias bandas del espectro
electromagnético. En esta oportunidad, presentamos nuestro proyecto para
estudiar tres galaxias vistas de canto, candidatas a presentar emisión por
fuera del plano galáctico en ondas de radio. Describiremos el criterio
utilizado para selecionar la muestra de galaxias y el software creado para
detectar la presencia de los vientos. Además, discutiremos los próximos pasos
a seguir de este proyecto.

Galaxies: halos Galaxies: starburst Radio continuum: galaxies Radiation
mechanisms: synchrotron

##  0.1 General

Starburst galaxies exhibit a remarkable level of star formation activity. The
combined effect of intense stellar winds generated by hot stars and a high
rate of supernova explosions leads to gas ejection from the star-forming
region into the surrounding environment, creating powerful winds (Chevalier &
Clegg, 1985). These winds play an important role in shaping the formation and
evolution of galaxies. They can inhibit star formation by depleting the
central gas reservoirs (Veilleux et al., 2005). As the gas escapes, star
formation decreases and the outflow stops. However, the arrival of fresh gas
can reignite the starburst. This influx of fresh gas is often triggered by
gravitational interactions with neighboring galaxies, or when previously
ejected gas that remains gravitationally bound falls back into the galaxy as
high-speed clouds (del Valle et al., 2018). However, some of the gas may not
be able to return to the galactic disk because it is no longer gravitationally
bound. Consequently, this gas enriches the intergalactic medium as it is
carried away by the outflow, supplementing it with heavy elements.

Different phases make up the outflow, and the relative positions of each phase
provide important clues to understanding the three-dimensional picture. The
hot gas in the cavity emits X-rays, the diffuse ionized gas surrounding the
hot gas emits $\rm H_{\alpha}$, and the atomic gas in a cooler shell pulled up
by a thick disk emits at the 21-cm line (see Fig. 8, Boosma et al., 2005). In
addition, the presence of magnetic fields and the relativistic particles that
emit synchrotron radiation can be studied by detecting the extra-planar radio
continuum emission (Krause et al., 2020). The best example is the starburst
galaxy NGC 253, located at a distance of 3.9 Mpc in the Sculptor group of
galaxies. The nuclear starburst is thought to be forming stars at a high rate
and producing a superwind that transports material into the halo. The galaxy
and its extra-planar gas have been well studied in recent years at various
wavelengths (Boosma et al., 2005; Bauer et al., 2007; Heesen et al., 2009a;
Romero et al., 2018).

We are carrying out a project to detect and characterize the extra-planar non-
thermal emission in three edge-on galaxies, to determine the non-thermal
contribution, to measure the magnetic field and its distribution, and to
detect embedded sources. In the following sections, we describe the sample of
galaxies, the method implemented to detect the continuum emission and the
current status of this project.

Table 1: Sample of galaxies Name | RA | DEC |  $i$ |  $D$ |  $\theta$ |  $\log(L_{\mathrm{IR}})$ |  $SFR$ |  $S_{60}/S_{100}$ |  $S_{1.4}$
---|---|---|---|---|---|---|---|---|---
| (hh:mm:ss) | (dd:mm:ss) | (deg) | (Mpc) | (arcmin) | ($L_{\sun}$ ) | ($\mathrm{M_{\sun}\leavevmode\nobreak\ yr^{-1}}$ ) |  | (mJy)
NGC 1055 | 02:41:45 | 00:26:35 | 83.8 | 17.7 | 6.17 | 10.09 | 1.83 | 0.36 | 200.9
NGC 4527 | 12:34:08 | 02:39:13 | 77.3 | 13.3 | 6.03 | 10.42 | 3.92 | 0.48 | 178
NGC 5690 | 14:37:41 | 02:17:27 | 83.3 | 19.14 | 3 | 10.23 | 2.53 | 0.4 | 72

##  0.2 The galaxy sample

From the Imperial IRAS-FSC Redshift Catalog (Wang & Rowan-Robinson, 2009), we
selected the galaxies according to the following criteria: 1. They should have
an inclination higher than $75^{\circ}$ to allow for a clear separation of the
halo emission from the disk emission. 2. They must fulfill that
$\log(L_{\mathrm{IR}}/L_{\sun})>10$ to guarantee the intense star formation
activity. 3. The ratio of infrared fluxes at 60 and 100 $\mu$m should be
$S_{60}/S_{100}>0.4$ to increase the chance of presenting disk-halo
interactions (Dahlem et al., 2001). 4. The angular sizes ($\theta$) and
distances ($D$) must satisfy $\theta>3$ arcmin and $D<20$ Mpc, to ensure their
detection at radio wavelengths and to allow us to resolve the emission in
z-direction.

The following galaxies constitute the galaxy sample. We list the main
properties of the galaxies, $i$ inclination, $D$ distance, $\theta$ angular
size, the $L_{\mathrm{IR}}$, SFR, $S_{60}/S_{100}$ and the flux density at 1.4
GHz ($S_{1.4}$) in Table LABEL:tab:sample.

NGC 1055: This galaxy shows evidence of a minor merger observed at the low-
angular resolution of the VLA and in Chandra X-ray observations (with low
integration time).

NGC 4527: This galaxy was detected in X-rays by Chandra. Low-angular
resolution observations carried out with the VLA suggest the presence of a
point-like source embedded in the diffuse radio continuum emission.

NGC 1421: This galaxy was previously observed with the VLA. Irwin et al.
(1999) suggest evidence of extra-planar emission and possible embedded
sources.

##  0.3 Methodology

Figure 1: Artificial image of a galaxy and its $z$-profile. Top: The
artificial image (upper panel) and the same image after rotation (lower
panel). The white contour level corresponds to the $3\sigma$ value. Bottom:
The $z$-profile for the artificial galaxy. The horizontal axis is the offset
from the galactic plane in arcseconds, while the vertical axis is the
logarithm of the flux density in $\mu$Jy/beam. The artificial data is shown in
magenta, the fitted disk is shown in dashed line and the halo component
corresponds to the dotted line. The total model is plotted in bold line.

We developed a software to identify and characterize the presence of extra-
planar emission. The software splits the galaxy image into multiple strips
perpendicular to the galactic disk (in $z$-direction) and averages them
together to create mean intensity profiles, hereafter $z$-profiles. The width
of this $z$-profile is given by the combined effect of the telescope
resolution, the galaxy inclination, the intrinsic height of the galactic
plane, and the possible presence of a radio halo. The software models these
components by assuming that the galactic plane is an infinitesimally thin
disk. When it is inclined and projected onto the sky, its profile along the
$z$-direction can be described either by a Gaussian or an exponential
function. Each function is convolved with a Gaussian-shaped beam to model the
effect introduced by the telescope resolution. If a single function is not
enough to fit the $z$-profile, a second exponential component is added to
model the emission from the halo.

We tested the software using artificial galaxies. These were created
considering different combinations of disk diameters and heights, besides an
inclination range between $80^{\circ}$ and $90^{\circ}$. Our software could
find the presence of the halos for artificial galaxies with inclinations
greater than $\sim 83^{\circ}$. We found no dependence of the halo scale
heights with the inclination. In Fig. 1, we show the results obtained for an
artificial galaxy created with an inclination of $84^{\circ}$ and a halo with
10 arcseconds FWHM. The best fit was found for a two-component function with
an exponential disk, properly identifying the halo.

##  0.4 Current status

Figure 2: NGC 4527 galaxy. The radio contour map at 550 MHz is superimposed on
a DSS image. The contour levels are 3.5, 5, 10, 20 in units of $\sigma$(=1 mJy
beam-1). The synthesised beam of $25^{\prime\prime}\times 25^{\prime\prime}$
is shown at the bottom-left corner.

We found that the software we developed to detect and characterize the extra-
planar emission works correctly by testing it using the artificial galaxies
described in Sect. 0.3. The three galaxies of the sample were observed with
the Giant Metrewave Radio Telescope (GMRT) at band-4 (550 MHz) and band-5
(1420 MHz). In order to attain 10$\mu$Jy beam-1 rms and considering a fudge
factor of 2.5 following ETC recommendations, we observed each source for 3 and
3.5 hours, at bands 5 and 4 respectively, including overhead. Because NGC 5690
is dimmer and the smallest, this galaxy was observed during 4 and 3.5 hours at
bands 5 and 4. Currently, we are calibrating and processing the data. A very
preliminary image of the galaxy NGC 4527 with synthesised beam of
$25^{\prime\prime}\times 25^{\prime\prime}$ and a noise level of 1 mJy beam-1,
is shown in Figure 2.

Acknowledgements: GER acknowledges financial support from the State Agency for
Research of the Spanish Ministry of Science and Innovation under grants
PID2019-105510GB-C31AEI/10.13039/501100011033/ and
PID2022-136828NB-C41/AEI/10.13039/501100011033/, and by“ERDF A way of making
Europ”, by the “European Unio”, and through the “Unit of Excellence María de
Maeztu 2020-2023” award to the Institute of Cosmos Sciences
(CEX2019-000918-M). Additional support came from PIP 0554 (CONICET).

## References

  * Bauer et al. (2007) Bauer, M., Pietsch, W., Trinchieri, G., et al. 2007, A&A, 467, 979
  * Boosma et al. (2005) Boosma, R., Oosterloo, T. A., Fraternali, F., et al. 2005, A&A, 431, 65
  * Chevalier & Clegg (1985) Chevalier, R. A. & Clegg, A. W. 1985, Nature, 317, 44
  * Dahlem et al. (2001) Dahlem, M., Lazendic, J. S., Haynes, R. F., Ehle, M., Lisenfeld, U. 2001, A&A, 374, 42
  * del Valle et al. (2018) del Valle, M. V., Müller, A. L., Romero, G. E. 2018, MNRAS, 475, 4298
  * Heesen et al. (2009a) Heesen, V., Beck, R., Krause, M., Dettmar, R. J. 2009a, A&A, 494, 563
  * Heesen et al. (2009b) Heesen, V., Krause, M., Beck, R., Dettmar, R. J. 2009b, A&A, 506, 1123
  * Irwin et al. (1999) Irwin, J. A., English, J., Sorathia, B. 1999, AJ, 117, 2102
  * Krause et al. (2020) Krause, M., Irwin, J., Schmidt, P., et al. 2020, A&A, 639, 112
  * Romero et al. (2018) Romero, G. E., Müller, A. L. & Roth, M. 2018, A&A, 616, A57
  * Veilleux et al. (2005) Veilleux , S., Cecil, G., Bland-Hawthorn, J. 2005, ARA&A, 43, 769
  * Wang & Rowan-Robinson (2009) Wang, L. & Rowan-Robinson, M. 2009, MNRAS, 398, 109


# Student $t$-Lévy regression model in YUIMA

Hiroki Masuda  Graduate School of Mathematical Sciences, University of Tokyo,
Japan<EMAIL_ADDRESS>,  Lorenzo Mercuri  Department of Economics,
Management and Quantitative Methods, University of Milan, Italy
INDAM-GNAMPA: Gruppo Nazionale per l’Analisi Matematica, la Probabilitá e le
loro Applicazioni
<EMAIL_ADDRESS>and  Yuma Uehara  Department of Mathematics, Faculty
of Engineering Science, Kansai University, Japan<EMAIL_ADDRESS>
###### Abstract.

The aim of this paper is to discuss an estimation and a simulation method in
the R package YUIMA for a linear regression model driven by a Student-$t$ Lévy
process with constant scale and arbitrary degrees of freedom. This process
finds applications in several fields, for example finance, physic, biology,
etc. The model presents two main issues. The first is related to the
simulation of a sample path at high-frequency level. Indeed, only the $t$-Lévy
increments defined on an unitary time interval are Student-$t$ distributed. In
YUIMA , we solve this problem by means of the inverse Fourier transform for
simulating the increments of a Student-$t$ Lévy defined on a interval with any
length. A second problem is due to the fact that joint estimation of trend,
scale, and degrees of freedom does not seem to have been investigated as yet.
In YUIMA , we develop a two-step estimation procedure that efficiently deals
with this issue. Numerical examples are given in order to explain methods and
classes used in the YUIMA package.

###### Key words and phrases:

Numerical implemetation, High-frequency sampling, Student Lévy Regression
Model.

##  1\. Introduction

YUIMA package in R provides several simulation and estimation methods for
stochastic processes [4, 16]. This paper introduces new classes and methods in
YUIMA for simulating and estimating a $t$-Lévy regression model based on high-
frequency observations. This model can be seen as a generalization of a
$t$-Lévy process (cf. [14, 7]) by adding covariates, which can be either
deterministic or stochastic processes. Adding covariates to a continuous-time
stochastic process is a widely used approach for constructing new processes in
various fields. For example, in finance, periodic deterministic covariates can
be employed to capture the seasonality observed in commodity markets [24].
Alternatively, in insurance and medicine, covariates are often incorporated
into mortality rate dynamics to account for age and cohort effects (see [13,
5] and references therein).

In such fields, data often exhibit heavy-tailed behavior in the margins and
thus, the inclusion of a $t$-Lévy process as a driving noise would be useful.
However, despite the simple mathematical definition of a $t$-Lévy regression
model, it poses significant challenges both in deriving its mathematical
properties and in implementing numerical methods. A key difficulty arises from
the fact that the $t$-Lévy driving noise is not closed under convolution.
Consequently, the distribution of its increments follows a $t$-distribution
only over a unit-time interval.

Motivated by this fact, in YUIMA , we propose three simulation methods for a
$t$-Lévy regression model. The key element is the construction of a random
number generator for the noise increments, which is essential for several
discretized simulation methods and involves the numerical inversion of the
characteristic function. More specifically, our method is based on the
approximation of the cumulative distribution function, and our method can
relieve numerical instability compared with directly using the density
function especially for simulating small time increments. Additionally, YUIMA
provides a two-stage estimation procedure and the corresponding estimator has
consistency and asymptotic normality as shown in [21].

The rest of this paper is organized as follows. We briefly summarize the
$t$-Lévy regression model in Section 2. We introduce the new YUIMA classes and
methods in Section 3 and we show some examples with simulated and real data
that highlight their usage in Section 4. Finally, Section 5 concludes the
paper.

##  2\. t-Lévy regression model

In this section, we review the main characteristics of the $t$-Lévy regression
model and the $t$-Lévy process used as its driving noise. We highlight the
main issues that arise in the simulation and estimation procedures for these
models. For a complete discussion on all mathematical aspects and the
properties of these procedures in both cases, we refer to [21].

The proposed $t$-Lévy regression model is a continuous-time stochastic process
of the form:

(2.1) $Y_{t}=X_{t}\cdot\mu+\sigma J_{t},\qquad t\in[0,T_{n}],$

where the $q$-dimensional vector process $X=(X_{t})$ with càdlàg paths
contains the covariates; the dot denotes the inner product in $\mathbb{R}^{q}$
and $J=(J_{t})$ is a $t$-Lévy process such that its unit-time distribution
$\mathcal{L}(J_{1})$ is given by:

(2.2) $\mathcal{L}(J_{1})=t_{\nu}:=t_{\nu}(0,1),$

where $t_{\nu}(\mu_{1},\sigma)$ denotes the scaled Student-$t$ distribution
with density:

$f(x;\mu,\sigma,\nu):=\frac{\Gamma(\frac{\nu+1}{2})}{\sigma\sqrt{\pi}\Gamma(\frac{\nu}{2})}\left\\{1+\left(\frac{x-\mu}{\sigma}\right)^{2}\right\\}^{-(\nu+1)/2}.$

The parameters $\nu>0$, $\mu_{1}\in\mathbb{R}$ and $\sigma>0$ represent the
degree of freedom, position, and scale parameters, respectively (see [14, 7]
for more details on a $t$-Lévy process). For this model, we consider the
situation where we estimate these three unknown parameters based a discrete-
time sample $\\{(X_{t_{j}},Y_{t_{j}})\\}_{j=0}^{[nT_{n}]}$ with
$t_{j}=t_{j}^{n}:=\frac{j}{n}$ and $T_{n}\to\infty$ as $n\to\infty$.

To develop the YUIMA simulation and estimation algorithms for the model in
(2.1), we need to address two problems: first, the simulation of the sample
path of $J=(J_{t})$ on a small time grid with $\Delta t\neq 1$ which is
essential for handling the increments:

$\Delta_{j}Y=\Delta_{j}X\cdot\mu+\sigma\Delta_{j}J,\quad j=0,1,\dots,[nT_{n}]$

where $\Delta_{j}Z$ denotes $Z_{t_{j}}-Z_{t_{j-1}}$ for any stochastic process
$Z=(Z_{t})$. Second, the identification of an efficient procedure for
estimating the model parameters $\left(\mu,\sigma\right)$ in (2.1) and the
degree of freedom $\nu$ in (2.2). We will introduce both the simulation and
estimation methods below.

Due to the stationary behavior of the Lévy increments, i.e.
$\mathcal{L}(\Delta_{i}J)=\mathcal{L}(J_{h})$ and
$h:=t_{i}-t_{i-1}=\frac{1}{n}$ for $i=1,\ldots,[nT_{n}]$, $\mathcal{L}(J_{h})$
admits the Lebesgue density:

$\displaystyle x$
$\displaystyle\mapsto\frac{1}{\pi}\int_{0}^{\infty}\cos(ux)\\{\varphi_{J_{1},\nu}(u)\\}^{h}du$
(2.3)
$\displaystyle=\left(\frac{2^{1-\nu/2}}{\Gamma(\nu/2)}\right)^{h}\frac{1}{\pi}\int_{0}^{\infty}\cos(ux)u^{\nu
h/2}\left(K_{\nu/2}(u)\right)^{h}du,$

where $\varphi_{J_{1},\nu}(u)$, the characteristic function of
$\mathcal{L}(J_{1})=t_{\nu}$ is given by

(2.4)
$\varphi_{J_{1},\nu}(u):=\frac{2^{1-\nu/2}}{\Gamma(\nu/2)}|u|^{\nu/2}K_{\nu/2}(|u|),\qquad
u\in\mathbb{R}$

and $K_{\nu}(t)$ denotes the modified Bessel function of the second kind
($\nu\in\mathbb{R}$, $t>0$):

$K_{\nu}(t)=\frac{1}{2}\int_{0}^{\infty}s^{\nu-1}\exp\left\\{-\frac{t}{2}\left(s+\frac{1}{s}\right)\right\\}ds.$

We here remark that although the explicit expression of (2.3) cannot be
obtained for general $h$, the tail index of $\mathcal{L}(J_{h})$ is the same
as that of $\mathcal{L}(J_{1})$ (cf. [3] Theorem 2), and hence, modeling tail
index via student $t$-Lévy process makes sense. A classical method for
simulating $t$-Lévy increments based on this Fourier inversion representation
is through the rejection method, as discussed in [15] and [8]. However, this
method can become time-consuming and unstable when dealing with very small
values of $h$, primarily due to the oscillatory behavior encountered during
the numerical evaluation of density. Such a problem often arises since to
express the high-frequently observed situation from a continuous process, we
should take a finer mesh of simulating the underlying process than that of
simulating the discrete observations. For this issue, in YUIMA , the Random
Number Generator for the increments is developed using the inverse quantile
method. The first step involves integrating the density in (2.3) to obtain the
cumulative distribution function (CDF). This approach leads to a more stable
behavior in the tails of the CDF compared to the density tails, owing to a
mitigation effect observed during the numerical evaluation of the following
double integral:

(2.5)
$F(y)=\int_{-\infty}^{y}\frac{1}{\pi}\int_{0}^{\infty}\cos(ux)\\{\varphi_{J_{1},\nu}(u)\\}^{h}dudx.$

Further details on our approach are provided in the next section. Another way
to simulate $t$-Lévy process is using series representation. Numerically, for
each time $t$, we need ad-hoc truncation of the infinite sum. The Gaussian
approximation of a small-jump part is also valid [1], but the associated error
may not be easy to control in a practical manner. For more details, see [19].

Following [21], YUIMA provides a two-stage estimation algorithm for
$\left(\mu,\sigma,\nu\right)\in\Theta_{\mu}\times\Theta_{\sigma}\times\Theta_{\nu}$.
Denote by $(\mu_{0},\sigma_{0},\nu_{0})$ the parameter true values of the
model in (2.1); the two-step algorithm can be summarized as follows:

  1. (1)

An estimator $\hat{a}:=(\hat{\mu}_{n},\hat{\sigma}_{n})$ of
$a_{0}=(\mu_{0},\sigma_{0})$ is obtained by maximizing the Cauchy quasi-
likelihood:

(2.6) $\hat{a}_{n}:=(\hat{\mu}_{n},\hat{\sigma}_{n})\in\mathop{\rm
argmax}_{a\in\overline{\Theta_{\mu}\times\Theta_{\sigma}}}\mathbb{H}_{1,n}(a).$

The Cauchy quasi-(log-)likelihood $\mathbb{H}_{1,n}(a)$ conditional on $X$ has
the following form:

$\displaystyle\mathbb{H}_{1,n}(a)$
$\displaystyle:=\sum_{j=1}^{N_{n}}\log\left\\{\frac{1}{h\sigma}\phi_{1}\left(\frac{\Delta_{j}Y-\mu\cdot\Delta_{j}X}{h\sigma}\right)\right\\}$
$\displaystyle=C_{n}-\sum_{j=1}^{N_{n}}\left\\{\log\sigma+\log\left(1+\epsilon_{j}(a)^{2}\right)\right\\},$

where $\phi_{1}$ denotes the density function of the standard Cauchy
distribution and the term $C_{n}$ does not depend on
$a=\left(\mu,\sigma\right)$. $N_{n}:=[nB_{n}]$ is the number of observations
in the part $[0,B_{n}]$ of the entire period $[0,T_{n}]$ where $(B_{n})$ is a
positive sequence satisfying $B_{n}\leq T_{n}$ and $\qquad
n^{\epsilon^{\prime\prime}}\lesssim B_{n}\lesssim n^{1-\epsilon^{\prime}}$ for
some $\epsilon^{\prime},\epsilon^{\prime\prime}\in(0,1)$.

  2. (2)

Then, we construct an estimator $\hat{\nu}_{n}$ of $\nu_{0}$ using the
Student-$t$ quasi-likelihood on the “unit-time” residual sequence
$\hat{\epsilon}_{i}$:

$\hat{\epsilon}_{i}:=\hat{\sigma}_{n}^{-1}\left(Y_{i}-Y_{i-1}-\hat{\mu}_{n}\cdot(X_{i}-X_{i-1})\right),$

for $i=1,\dots,[T_{n}]$, which is expected to be approximately i.i.d.
$t_{\nu}$-distributed. Therefore $\hat{\nu}_{n}$ solves:

$\hat{\nu}_{n}\in\mathop{\rm
argmax}_{\nu\in\overline{\Theta}_{\nu}}\mathbb{H}_{2,n}(\nu),$

where

$\mathbb{H}_{2,n}(\nu):=\sum_{i=1}^{[T_{n}]}\left(-\frac{1}{2}\log\pi+\log\Gamma\left(\frac{\nu+1}{2}\right)-\log\Gamma\left(\frac{\nu}{2}\right)-\frac{\nu+1}{2}\log\left(1+\hat{\epsilon}_{i}^{2}\right)\right).$

We note that the first estimation scheme is based on the locally Cauchy
property of $J$: $h^{-1}J_{h}\xrightarrow{\mathcal{L}}t_{1}$ (standard Cauchy)
as $h\to 0$. Let $\hat{u}_{a,n}:=\sqrt{N_{n}}(\hat{a}_{n}-a_{0})$ and
$\hat{u}_{\nu,n}:=\sqrt{T_{n}}(\hat{\nu}_{n}-\nu_{0})$. Under some regularity
conditions on the covariate process $X=\left(X_{t}\right)$ (see Assumption 2.1
in [21]for all requirements on $X$), the estimators have the following joint
asymptotic normality:

$(\hat{\Gamma}_{a,n}^{1/2}\hat{u}_{a,n},\hat{\Gamma}_{\nu,n}^{1/2}\hat{u}_{\nu,n})\xrightarrow{\mathcal{L}}N_{q+2}(0,I_{q+2}),$

where $\psi_{1}$ denotes the trigamma function and

$\displaystyle\hat{\Gamma}_{a,n}={\rm
diag}\left(\frac{1}{2\hat{\sigma}_{n}^{2}N_{n}}\sum_{j=1}^{N_{n}}\left(\frac{1}{h}\Delta_{j}X\right)^{\otimes
2},\frac{1}{2\hat{\sigma}_{n}^{2}}\right),$
$\displaystyle\hat{\Gamma}_{\nu,n}=\frac{1}{4}\Bigg{(}\psi_{1}\Bigg{(}\frac{\hat{\nu}_{n}}{2}\Bigg{)}-\psi_{1}\Bigg{(}\frac{\hat{\nu}_{n}+1}{2}\Bigg{)}\Bigg{)}.$

Since $\hat{\Gamma}_{a,n}^{1/2}$ and $\hat{\Gamma}_{\nu,n}$ can be constructed
only by the observations, we can easily obtain the confidence intervals of
each parameter.

##  3\. Classes and Methods for t-Lévy Regression Models

This section provides an overview of the new classes and methods introduced in
YUIMA for the mathematical definition, trajectory simulation, and estimation
of a Student Lévy Regression model. To handle this model, the first step
involves constructing an object of the yuima.LevyRM-class. As an extension of
the yuima-class (refer to [4] for more details), yuima.LevyRM-class inherits
slots such as @data, @model, @sampling, and @functional from its parent class.
The remaining slots store specific information related to the Student-$t$-Lévy
regression model.

Notably, the slot `@unit_Levy` contains an object of the yuima.th class, which
represents the mathematical description of the Student-$t$ Lévy process
$J_{t}$ (see the subsequent section for detailed explanations). The labels of
the regressors are saved in the slot @regressors, while slots @LevyRM and
@paramRM respectively cache the names of the output process $Y_{t}$ and a
string vector reporting the regressors’ coefficients, the scale parameter, and
the degree of freedom. The yuima.th-class is obtained by the new `setLaw_th`
constructor. This function requires the arguments used for the numerical
inversion of the characteristic function, and its usage is discussed in the
next section.

Once the yuima.th-object is created, we define the system of stochastic
differential equations (SDEs) that describes the behavior of the regressors,
with their mathematical definitions stored in an object of the yuima.model-
class. Both yuima.th and yuima.model objects are used as inputs for the setLRM
constructor, which returns an object of the yuima.LevyRM-class. The following
chunk code reports the input for this new function.

setLRM(unit_Levy, yuima_regressors, LevyRM = "Y", coeff = c("mu", "sigma0"),
data = NULL, sampling = NULL, characteristic = NULL, functional = NULL) As is
customary for any class extending the yuima-class, the simulate method enables
the generation of sample paths for the Student Lévy Regression model. To
simulate trajectories, an object of the yuima.sampling-class is constructed to
represent an equally spaced grid-time used in trajectory simulation. The
regressors’ paths are obtained using the Euler scheme, while the increments of
the Student-$t$ Lévy process are simulated using the random number generator
available in the slot @rng of the yuima.th-object.

The last method available in YUIMA is `estimation_LRM`. This method allows the
users to estimate the model using either real or simulated data. The
estimation follows a two-step procedure, introduced in [21].

estimation_LRM(start, model, data, upper, lower) For this function, the
minimal inputs are start, model, data, upper and lower. The arguments start
are the initial points for the optimization routine while upper and lower
corresponds to the box constraints. The yuima.LevyRM-object is passed to the
function through the input model while the input data is used to pass the
dataset to the internal optimization routine. Figure 1 describes these new
classes and methods, along with their respective usage.

yuima.LevyRM setLRM yuima.th yuima.model `setLaw_th` Inputs for the Fourier
Inv. simulate yuima.sampling Sample Path `estimation_LRM` yuima.data Estimated
Param. Figure 1. Scheme of classes and methods in YUIMA for the t-Lévy
Regression model.

###  3.1. yuima.th: A new class for mathematical description of a Student-t
Lévy process

In this section, the steps for the construction of an yuima.th-object are
presented. As remarked in Section 3, this object contains all information on
the Student-$t$ Lévy process. Moreover, detailed information is provided
regarding the numerical algorithms utilized for evaluating the density
function (2.3). The construction of this object is accomplished using the
`setLaw_th` constructor, and the subsequent code snippet displays its
corresponding arguments.

setLaw_th(h = 1, method = "LAG", up = 7, low = -7, N = 180, N_grid = 1000,
regular_par = NULL) The input h is the length of the step-size of each time
interval for the Student-$t$ Lévy increment $\Delta J_{h}=J_{h}-J_{0}$. Its
default value, h=1, indicates that the yuima.th-object describes completely
the process $J_{t}$ at time 1. The argument method refers to the type of
quadrature used for the computation of the integral in (2.3) while the
remaining arguments govern the precision of the integration routine.

The yuima.th-class inherits the slots @rng, @density, @cdf and @quantile from
its parent class yuima-law [20]. These slots store respectively the random
number generator, the density function, the cumulative distribution function
and the quantile function of $\Delta J_{h}$.
As mentioned in Section 1, the density function with $h\neq 1$ does not have a
closed-form formula and, therefore, the inversion of the characteristic
function is necessary. YUIMA provides three methods for this purpose: the
Laguerre quadrature, the COS method and the Fast Fourier Transform. The Gauss-
Laguerre quadrature is a numerical integration method employed for evaluation
integrals in the following form:

$\mathcal{I}=\int_{0}^{+\infty}f\left(x\right)e^{-x}\mbox{d}x.$

This procedure has been recently used for the computation of the density of
the variance gamma and the transition density of a CARMA model driven by a
time-changed Brownian motion [18, 22], this motivates its application in this
paper.

Let $f\left(x\right)$ be a continuous function defined on
$\left[0,+\infty\right)$ such that:

$\mathcal{I}=\int_{0}^{+\infty}f\left(x\right)e^{-x}\mbox{d}x<+\infty;$

the integral $\mathcal{I}$ can be approximated as follows:

(3.1)
$\mathcal{I}\approx\sum_{j=1}^{N}\omega\left(k_{j}\right)f\left(k_{j}\right),$

where $k_{j}$ is the $j$th-root of the $N$-order Laguerre polynomial111The
Laguerre polynomial can be defined recursively as follows:
$L_{0}\left(x\right),L_{1}(x)=1-x,\ldots,L_{N}\left(x\right)=\frac{\left(2N-x\right)L_{N-1}\left(x\right)-\left(N-1\right)L_{N-2}\left(x\right)}{N}$
$L_{N}\left(x\right)$ and the weights $\omega\left(k_{j}\right),j=1,\ldots,N$
are defined as:

(3.2)
$\omega\left(k_{j}\right)=\frac{k_{j}}{\left(N+1\right)^{2}L_{N+1}^{2}\left(k_{j}\right)}.$

To apply the approximation in (3.1), we rewrite the inversion formula in (2.3)
as follows:

(3.3) $\displaystyle f\left(x\right)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\left(\frac{2^{1-\frac{\nu}{2}}}{\Gamma\left(\nu/2\right)}\right)\int_{0}^{+\infty}\cos\left(ux\right)u^{\nu
h/2}\left(K_{\nu/2}\left(u\right)\right)^{h}\mbox{d}u$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\left(\frac{2^{1-\frac{\nu}{2}}}{\Gamma\left(\nu/2\right)}\right)\int_{0}^{+\infty}\cos\left(ux\right)u^{\nu
h/2}\left(K_{\nu/2}\left(u\right)\right)^{h}e^{u}e^{-u}\mbox{d}u.$

Applying the result in (3.1), the density function of $J_{h}$ can be
approximated with the formula reported below:

$\hat{f}_{N}\left(x_{j}\right)=\frac{1}{\pi}\left(\frac{2^{1-\frac{\nu}{2}}}{\Gamma\left(\nu/2\right)}\right)\sum_{j=1}^{N}\cos\left(k_{j}x\right)k_{j}^{\nu
h/2}\left(K_{\nu/2}\left(k_{j}\right)\right)^{h}e^{k_{j}}\omega\left(k_{j}\right).$

Notably, the approximation formula’s precision in equation (3.1) can be
enhanced through the argument N in the `setLaw_th` constructor. The roots
$k_{j}$ and the weights $\omega\left(k_{j}\right)$ are internally computed
using the gauss.quad function from the R package statmod, with a maximum
allowed order of 180 for the Laguerre polynomial.

The COS method is based on the Fourier Cosine expansion employed for an even
function with the compact domain $\left[-\pi,\pi\right]$. This method has been
widely applied in the finance literature for the computation of the exercise
probability of an option and its no-arbitrage price, we refer to [10, 11] and
references therein for details. Let
$g\left(\theta\right):\left[-\pi,\pi\right]\rightarrow\mathbb{R}$ be an even
function, its Fourier Cosine expansion reads:

(3.4)
$g\left(\theta\right)=\frac{1}{2}a_{0}+\sum_{k=1}^{\infty}a_{k}\cos\left(k\theta\right)$

with

$a_{k}=\frac{2}{\pi}\int_{0}^{\pi}g\left(\theta\right)\cos\left(k\theta\right)\mbox{d}\theta.$

Denoting with $f\left(x\right)$ the density function of the increment $\Delta
J_{h}$, the new function $\bar{g}\left(\theta\right)$ is defined as:

(3.5)
$\bar{g}\left(\theta\right):=f\left(\frac{L}{\pi}\theta\right)\frac{L}{\pi}\mathbbm{1}_{\left\\{-\pi\leq\theta\leq\pi\right\\}},$

The function $\bar{g}\left(\theta\right)$ is still an even function for any
$L>0$ and the Fourier Cosine expansion in (3.4) can be applied. The
coefficient $a_{k}$ is determined as follows:

$\displaystyle a_{k}$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\int_{0}^{\pi}\bar{g}\left(\theta\right)\cos\left(k\theta\right)\mbox{d}\theta$
$\displaystyle=$
$\displaystyle\frac{2}{\pi}\int_{0}^{\pi}f\left(\frac{L}{\pi}\theta\right)\cos\left(k\theta\right)\frac{L}{\pi}\mbox{d}\theta.$

Setting $\theta=\frac{\pi}{L}x$, we have:

(3.6)
$a_{k}=\frac{2}{\pi}\int_{0}^{L}f\left(x\right)\cos\left(k\frac{\pi}{L}x\right)\mbox{d}x.$

The coefficient $a_{k}$ can be rewritten using the characteristic function of
$\Delta J_{h}$ at $k\frac{\pi}{L}$, i.e.:

(3.7)
$a_{k}=\frac{2}{\pi}\left[\varphi\left(k\frac{\pi}{L}\right)-\int_{L}^{+\infty}f\left(x\right)\cos\left(k\frac{\pi}{L}x\right)\mbox{d}x\right].$

For a sufficiently large $L$, the coefficient $a_{k}$ can be approximated as
follows:

(3.8) $a_{k}\approx\frac{2}{\pi}\varphi\left(k\frac{\pi}{L}\right).$

The following series expansion is achieved:

(3.9)
$f\left(x\right)=\frac{1}{2}\bar{a}_{0}+\sum_{k=1}^{+\infty}\bar{a}_{k}\cos\left(k\frac{\pi}{L}x\right)$

with

$\bar{a}_{k}=a_{k}\frac{\pi}{L}\approx\frac{2}{L}\varphi\left(k\frac{\pi}{L}\right).$

Finally, (3.9) can be approximated by truncating the series as follows:

(3.10)
$\hat{f}_{N}\left(x\right)=\frac{1}{2}\hat{a}_{0,L}+\sum_{k=1}^{N}\hat{a}_{k,L}\cos\left(k\frac{\pi}{L}x\right)$

with

$\hat{a}_{k,L}=\frac{2}{L}\varphi\left(k\frac{\pi}{L}\right).$

The precision of $\hat{f}_{N,L}\left(x\right)$ in (3.10) depends on $N$ and
$L$. Users can select these two quantities using N and the couple (up, low)
respectively. $L$ is computed internally in the `setLaw_th` constructor as L =
max(|low|,up).

The last method available in YUIMA is the Fast Fourier Transform (FFT) [23, 6]
for the inversion of the characteristic function. YUIMA uses the FFT method
developed internally in the function `FromCF2yuima_law` for the inversion of
any characteristic function defined by users. In the latter case, the density
function of $\Delta J_{h}$ is based on the following general inversion
formula:

(3.11)
$f\left(x\right)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{-iux}\varphi\left(u\right)\mbox{d}u.$

As first step, we apply to the integral in (3.11) the change variable
$u=2\pi\omega$ and we have:

(3.12) $f\left(x\right)=\int_{-\infty}^{+\infty}e^{-i2\pi\omega
x}\varphi\left(2\pi\omega\right)\mbox{d}\omega.$

We consider discrete support for $x$ and $\omega$. The $x$ grid has the
following structure:

(3.13) $x_{0}=a,\ldots,x_{j}=a+j\Delta x,\ldots,x_{N}=b$

with $\Delta x=\frac{b-a}{N}$ where $N+1$ is the number of the points in the
grid in (3.13). Similarly, we define a $\omega$ grid:

(3.14)
$\omega_{0}=-\frac{N}{2}\Delta\omega,\ldots,\omega_{n}=-\frac{N}{2}+n\Delta\omega\ldots,\omega_{N}=\frac{N}{2}\Delta\omega$

where $\Delta\omega=\frac{1}{N\Delta x}=\frac{1}{b-a}$. Both grids have the
same dimension $N+1$, $\Delta x$ shrinks as $N\rightarrow+\infty$ while
$\Delta\omega$ reduces for large values of $b-a$. For any $x_{j}$ in the grid
we can approximate the integral in (3.11) using the left Riemann summation,
therefore, we have the approximation $\hat{f}_{N}\left(x_{j}\right)$:

$\displaystyle\hat{f}_{N}\left(x_{j}\right)$
$\displaystyle=\Delta\omega\sum_{n=1}^{N}e^{2i\pi\left(-\frac{N}{2}\Delta\omega+\left(n-1\right)\Delta
u\right)x_{j}}\varphi\left(-\pi N\Delta u+2\pi\left(n-1\right)\right)$
$\displaystyle=\Delta\omega e^{i\pi
N\Delta\omega}\sum_{n=1}^{N}e^{-2i\pi\left(n-1\right)\Delta
u\left(a+\left(j-1\right)\Delta x\right)}\varphi\left(-\pi
N\Delta\omega+2\pi\left(n-1\right)\right)$ (3.15) $\displaystyle=\Delta\omega
e^{i\pi
N\Delta\omega}\sum_{n=1}^{N}e^{\frac{-2i\pi\left(n-1\right)\left(j-1\right)}{N}}\varphi\left(-\pi
N\Delta\omega+2\pi\left(n-1\right)\right)e^{\frac{-2i\pi\left(n-1\right)a}{N}}.$

The last equality in (3.15) is due to the identity $\Delta\omega\Delta
x=\frac{1}{N}$. To evaluate the summation in (3.15), we use the FFT algorithm
and

(3.16) $\hat{f}_{N}\left(x_{j}\right)=\Delta\omega e^{i\pi
N\Delta\omega}\texttt{FFT}\left[\varphi\left(-\pi
N\Delta\omega+2\pi\left(n-1\right)\right)e^{\frac{-2i\pi\left(n-1\right)a}{N}}\right].$

In this case, we have two sources of approximation errors that we can control
using the arguments up, low and N in the `setLaw_th` constructor. N denotes
the number of intervals in the grid used for the Left-Riemann summation while
up, low are used to compute the step size of this grid.

Once the density function has been obtained using one of the three methods
described above, it is possible to approximate the cumulative distribution
function using the Left-Riemann summation computed on the grid in (3.13),
therefore the cumulative distribution function $F\left(\cdot\right)$ for each
$x_{j}$ on the grid is determined as follows:

(3.17)
$\hat{F}\left(x_{j}\right)=\sum_{x_{k}<x_{j}}\hat{f}_{N}\left(x_{k}\right)\Delta
x.$

In this way, we can construct a table that we can use internally in the cdf
and quantile functions. Moreover, to evaluate the cumulative distribution
function at any $x\in\left(x_{j-1},x_{j}\right)$, we interpolate linearly its
value using the couples $\left(x_{j-1},\hat{F}\left(x_{j-1}\right)\right)$ and
$\left(x_{j},\hat{F}\left(x_{j}\right)\right)$. The random numbers can be
obtained using the inversion sampling method.

We conclude this section with a numerical comparison among the three methods
implemented in YUIMA . To assess the precision of our code we use as a target
the cumulative distribution function of a Student-$t$ with $\nu=3$ computed
through the R function pt available in the package stats. To conduct this
comparison, we construct three yuima.th-objects as displayed in the code
snippet reported below.

# To instal YUIMA from Github repositoryR> library(devtools)R>
install_github("yuimaproject/yuima")########### Inputs ###########R>
library(yuima)R> nu <\- 3R> h <\- 1 # step size for the interval timeR> up <\-
10R> low <\- -10# Support definition for variable xR> x <\-
seq(low,up,length.out=100001)# Definition of yuima.th-objectR> law_LAG <\-
setLaw_th(h = h, method = "LAG", up = up, low = low, N = 180) # LaguerreR>
law_COS <\- setLaw_th(h = h, method = "COS", up = up, low = low, N = 180) #
COSR> law_FFT <\- setLaw_th(h = h, method = "FFT", up = up, low = low, N =
180) # FFT# Cumulative Distribution Function: we apply the cdf method to the
yuima.th-objectR> Lag_time <\- system.time(ycdf_LAG <\- cdf(law_LAG, x,
list(nu=nu))) # LaguerreR> COS_time <\- system.time(ycdf_COS<-cdf(law_COS, x,
list(nu=nu))) # COSR> FFT_time <\- system.time(ycdf_FFT<-cdf(law_FFT, x,
list(nu=nu))) # FFT User can specify the numerical methods of the inversion of
the characteristic function using the input method. This argument assumes
three values: "LAG" for the Gauss-Laguerre quadrature, "COS" for the Cosine
Series Expansion and "FFT" for the Fast Fourier Transform. Once the yuima.th-
object has been constructed, its cumulative distribution function is computed
by applying the YUIMA method222An object of yuima.th-class inherits the dens
method for the density computation, cdf for the evaluation of the cumulative
distribution function, quantile for the quantile function and rand for the
generation of the random sample. We refer to [20] for the usage of these
methods. cdf that returns the cumulative distribution function for the numeric
vector x.

Figure 2. Cumulative distribution function comparison for a Student-$t$ random
variable with $\nu=3$.

Figure 2 compares the cumulative distribution functions of the Student-$t$
obtained using the methods in YUIMA with the pt function. In all cases, we
have a good level of precision that is also confirmed by Table 1. As expected
the fastest method is the FFT which seems to be also the most precise.

| COS | FFT | LAG
---|---|---|---
RMSE | 0.021 | 0.021 | 0.064
Max | 0.032 | 0.032 | 0.085
Min | 7.88e-05 | 2.18e-05 | 4.97e-04
sec. | 1.47 | 4.00e-02 | 1.2
Table 1. Summary comparison of the distance between each method available in
YUIMA and the cumulative distribution function obtained using the R function
pt.

To further investigate this fact, we study the behaviour of the cumulative
distribution function when $h$ varies. For $h=0.01$, we observe an oscillatory
behaviour on tails for the FFT and COS while the Laguerre method seems to be
more stable as shown in Figure 3. To have a fair comparison, we set N = 180,
however we notice that the precision of FFT can be drastically improved by
tuning the inputs N, up and low.

Figure 3. Cumulative distribution function comparison for a Student-$t$ Lévy
increment with $\nu=3$ and $h=0.01$.

##  4\. Numerical examples

This section presents a series of numerical examples demonstrating the
practical application of newly developed classes and methods within the
Student-$t$ Regression model. Specifically, we showcase the simulation and
estimation of models where the regressors are determined by deterministic
functions of time in the first example. The second example introduces
integrated stochastic regressors. Lastly, we perform an analysis using real
data in the final example.

###  4.1. Model with deterministic regressors.

In this example, we consider two deterministic regressors and the dynamics of
the model have the following form:

(4.1) $Y_{t}=\mu_{1}\cos\left(5t\right)+\mu_{2}\sin\left(t\right)+\sigma
J_{t}$

with the true values
$\left(\mu_{1},\mu_{2},\sigma_{0},\nu_{0}\right)=\left(5,-1,3,3\right)$. The
estimation of the model (4.1) has been investigated in [21] where the
empirical distribution of the studentized estimates has been discussed. To use
the simulation method in YUIMA we have to write the dynamics of the regressors
$X_{1,t}:=\cos\left(5t\right)$ and $X_{2,t}=:\sin\left(5t\right)$, that is:

(4.2) $\mbox{d}\left[\begin{array}[]{c}X_{1,t}\\\
X_{2,t}\end{array}\right]=\left[\begin{array}[]{c}-5\sin\left(5t\right)\\\
\cos\left(t\right)\end{array}\right]\mbox{d}t$

with the initial condition:

$X_{1,0}=1,\ X_{2,0}=0.$

In the following, we show how to implement this model in YUIMA . In this
example, we set $h_{n}=1/50$, and thus the number of the observations $n$ over
unit time is $50$. We also set the terminal time of the whole observations
$T_{n}=50$ and that of the partial observations $B_{n}=15$.

R> library(yuima)########### Inputs ############ Inputs for Fourier
InversionR> method_Fourier = "FFT"; up = 6; low = -6; N = 2^17; N_grid =
60000# Inputs for the sample grid timeR> Factor <\- 1R> Factor1 <\- 1R>
initial <\- 0; Final_Time <\- 50 * Factor; h <\- 0.02/Factor1

We notice that the variable Factor1 controls the step size of the time grid.
Indeed for different values of this quantity, we have a different value for
$h$ for example Factor1 = 1, 2, 4, 7.2 corresponds h = 0.02, 0.01, 0.005,
1/365.

The first step is the definition of an object of yuima.th-class using the
constructor `setLaw_th`. This object contains the random number generator, the
density function, the cumulative distribution function and the quantile
function for constructing the increments of a Student-$t$ Lévy process.

######################################## Example 1: Deterministic Regressors
########################################R> mu1 <\- 5; mu2 <\- -1; scale <\- 3;
nu <\- 3 # Model Parameters# Model DefinitionR> law1 <\- setLaw_th(method =
method_Fourier, up = up, low = low, N = N, N_grid = N_grid) # yuima.thR>
class(yuima_law)

[1] "yuima.th"attr(,"package")[1] "yuima"

The next step is to define the dynamics of the regressors described in (4.2).
This set of differential equations is defined in YUIMA using the standard
constructor setModel. Once an object containing the mathematical description
of the regressors has been defined, we use setLRM to obtain the yuima.LevyRM.
In the following, we report the command lines for the definition of the
Student Lévy Regression Model in (4.1).

R> regr1 <\- setModel(drift = c("-5*sin(5*t)", "cos(t)"), diffusion =
matrix("0",2,1), solve.variable = c("X1","X2"), xinit = c(1,0)) # Regressors
definitionR> Mod1 <\- setLRM(unit_Levy = law1, yuima_regressors = regr1) #
t-Regression ModelR> class(Mod1)

[1] "yuima.LevyRM"attr(,"package")[1] "yuima"

Using the object Mod1 we simulate a trajectory of the model in (4.1) using the
YUIMA method simulate.

# SimulationR> samp <\- setSampling(initial, Final_Time, n = Final_Time/h)R>
true.par <\- unlist(list(mu1 = mu1, mu2 = mu2, sigma0 = scale, nu = nu))R>
set.seed(1)R> sim1 <\- simulate(Mod1, true.parameter = true.par, sampling =
samp) Figure 4 reports the simulated sample paths for the regressors
$X_{1,t}$, $X_{2,t}$ and the Student Lévy Regression model $Y_{t}$.

Figure 4. Simulated Trajectory of the Student Lévy Regression model defined in
(4.1).

The next step is to study the behaviour of the two-step estimation procedure
described in Section 2. To run this procedure, we use the method
`estimation_LRM` and then we initialize randomly the optimization routine as
shown in the following command lines.

# EstimationR> lower <\- list(mu1 = -10, mu2 = -10, sigma0 = 0.1)R> upper <\-
list(mu1 = 10, mu2 = 10, sigma0 = 10.01)R> start <\- list(mu1 = runif(1, -10,
10), mu2 = runif(1, -10, 10),\+ sigma0 = runif(1, 0.01, 4))R> Bn <\-
15*FactorR> est1 <\- estimation_LRM(start = start, model = sim1, data =
sim1@data,\+ upper = upper, lower = lower, PT = Bn)R> class(est1)

[1] "yuima.qmle"attr(,"package")[1] "yuima"

R> summary(est1)

Quasi-Maximum likelihood estimationCall:estimation_LRM(start = start, model =
sim1, data = sim1@data, upper = upper, lower = lower, PT = Bn)Coefficients:
Estimate Std. Errormu1 5.029555 0.03538355mu2 -1.128905 0.18026088sigma0
2.425147 0.12523404nu 2.735344 0.47457435-2 log L: 3236.535 73.3833

It is also possible to construct the dataset without the method simulate. This
result can be achieved by constructing an object of yuima.th-class that
represents the increment of the Student-$t$ Lévy process on the time interval
with length $h$.

# Dataset construction without YUIMAR> time_t <\-
index(get.zoo.data(sim1@data)[[1]])R> X1 <\- zoo(cos(5*time_t), order.by =
time_t)R> X2 <\- zoo(sin(time_t), order.by = time_t)R> law1_h <\- setLaw_th(h
= h, method = method_Fourier, up = up, low = low, N = N, N_grid = N_grid)R>
print(c(law1_h@h,law1@h))

[1] 0.02 1.00 The object `law1_h` refers to the Student-$t$ Lévy increment
over an interval with length $0.02$ as it can be seen looking at the slot
...@h.

R> set.seed(1)R> names(nu)<\- "nu"R> J_t <\- zoo(cumsum(c(0,rand(law1_h,
Final_Time/h,nu))), order.by = time_t)R> Y <\- mu1 * X1 + mu2 * X2 +
scale/sqrt(nu) * J_tR> data1_a <\- merge(X1, X2, Y)

The estimation is performed by means of the method `estimation_LRM` as in the
previous example.

R> est1_a <\- estimation_LRM(start = start, model = Mod1, data =
setData(data1_a),\+ upper = upper, lower = lower, PT = Bn)R> summary(est1_a)

Quasi-Maximum likelihood estimationCall:estimation_LRM(start = start, model =
Mod1, data = setData(data1_a), upper = upper, lower = lower, PT =
Bn)Coefficients: Estimate Std. Errormu1 4.987778 0.03642823mu2 -1.164138
0.18558422sigma0 2.495930 0.12888927nu 2.407683 0.41221590-2 log L: 3232.784
85.32438

Using the result stored in the summary we can construct a confidence interval
using the command lines reported below.

R> info_sum <\- summary(est1_a)@coefR> alpha <\- 0.025R> Confidence_Int_95 <\-
rbind(info_sum[,1]+info_sum[,2]*qnorm(alpha),\+
info_sum[,1]+info_sum[,2]*qnorm(1-alpha),unlist(true.par),\+ info_sum[,1])R>
rownames(Confidence_Int_95) <\- c("LB","UB","True_par","Est_par")R>
print(Confidence_Int_95, digit = 4)

mu1 mu2 sigma0 nuLB 4.916 -1.5279 2.243 1.600UB 5.059 -0.8004 2.749
3.216True_par 5.000 -1.0000 3.000 3.000Est_par 4.988 -1.1641 2.496 2.408 Based
on the aforementioned examples, the estimated value for the $\sigma$ parameter
deviates from its true value. To further investigate this observation, we
conducted a comparative analysis using the three integration methods discussed
in Section 3.1. This exercise was repeated for three different values of
$h=0.01,0.005$ and $\frac{1}{365}$. To ensure a fair comparison among the
Laguerre, Cosine Series expansion, and Fast Fourier Transform methods, we set
the argument N = 180. The obtained results are presented in Table 2.

| Laguerre Method | COS Method | FFT Method
---|---|---|---
$h$ | 0.01 | 0.005 | $1/365$ | 0.01 | 0.005 | $1/365$ | 0.01 | 0.005 | $1/365$
$\hat{\mu_{1}}$ | 4.968 | 5.019 | 4.976 | 4.939 | 5.065 | 4.875 | 4.934 | 5.069 | 4.876
| (0.029) | (0.020) | (0.014) | (0.041) | (0.049) | (0.054) | (0.042) | (0.049) | (0.054)
$\hat{\mu_{2}}$ | -1.065 | -1.044 | -0.914 | -1.192 | -1.195 | -0.555 | -1.179 | -1.178 | -0.528
| (0.146) | (0.104) | (0.074) | (0.210) | (0.247) | (0.277) | (0.213) | (0.248) | (0.277)
$\hat{\sigma_{0}}$ | 2.770 | 2.794 | 2.657 | 3.993 | 6.654 | 10.010 | 4.055 | 6.674 | 10.010
| (0.101) | (0.072) | (0.051) | (0.145) | (0.172) | (0.192) | (0.148) | (0.172) | (0.193)
$\hat{\nu}$ | 3.462 | 3.281 | 3.067 | 2.827 | 2.223 | 2.227 | 5.749 | 8.310 | 15.175
| (0.615) | (0.580) | (0.538) | (0.492) | (0.377) | (0.378) | (1.063) | (1.571) | (2.940)
sec. | 2.07 | 2.02 | 2.07 | 1.10 | 1.24 | 1.34 | 0.15 | 0.24 | 0.39
Table 2. Estimated parameters for $h=0.01,0.005,1/365$ and different
integration methods. The number of points for the inversion of the
characteristic function is $N=180$. The parenthesis shows the asymptotic
standard error. The last row reports the seconds necessary for the simulation
of a sample path.

In the Laguerre method, we observe that reducing the step size $h$ leads to
improved estimates of $\nu$, however looking at the asymptotic standard error,
the estimates for $\sigma_{0}$ seem to maintain the bias. As for the remaining
methods, the estimates exhibit unsatisfactory performance due to the imposed
restriction of N = 180 (that is the maximum value allowed for the Laguerre
quadrature in the evaluation of the density in (3.3)). This outcome is not
unexpected, given the Laguerre method’s ability to yield a valid distribution
even with a relatively small value of N, as demonstrated in Figure 3.
Nevertheless for the COS and the FFT, it is possible to enhance the results by
increasing the value of the argument N that yields a more accurate result in
the simulation of the sample path for the model described in (4.1).

| COS Method | FFT Method
---|---|---
N | 1000 | 5000 | 10000 | 1000 | 5000 | 10000
$\mu_{1}$ | 4.968 | 4.975 | 4.975 | 4.968 | 4.975 | 4.975
| (0.034) | (0.016) | (0.016) | (0.035) | (0.306) | (0.016)
$\mu_{2}$ | -0.887 | -0.909 | -0.909 | -0.886 | -0.908 | -0.908
| (0.174) | (0.082) | (0.082) | (0.177) | (0.022) | (0.082)
$\sigma_{0}$ | 3.300 | 2.964 | 2.964 | 3.363 | 2.963 | 2.964
| (0.064) | (0.057) | (0.057) | (0.065) | (0.057) | (0.057)
$\nu$ | 2.713 | 2.765 | 2.765 | 3.279 | 2.760 | 2.702
| (0.470) | (0.480) | (0.480) | (0.579) | (0.479) | (0.468)
sec. | 3.49 | 14.83 | 33.70 | 0.38 | 0.44 | 0.45
Table 3. Estimated parameters for $N=1000,5000,10^{4}$ in the COS and the FFT
methods. The step size $h$ is $1/365$. The parenthesis shows the asymptotic
standard error. The last row reports the seconds necessary for the simulation
of a sample path.

Table 3 shows the estimated parameters for the varying value of N and
$h=1/365$. As expected increasing the precision in the quadrature improved
estimates for both methods. Notably for N $\geq$ 5000, all estimates fall
within the asymptotic confidence interval at the $95\%$ level.

###  4.2. Model with integrated stochastic regressors

In this section, we consider two examples whose regressors are stochastic.
Moreover, to satisfy the regularity conditions, the regressors are supposed to
be an integrated version of stochastic processes.

#### Example 1

In this example, we consider the following continuous time regression model
with a single regressor:

(4.3) $Y_{t}=\mu\int_{0}^{t}X_{u}du+\sigma J_{t},\quad J_{1}\sim t_{\nu},$

with the true values $(\mu_{0},\sigma_{0},\nu_{0})=(-3,3,2.5)$, and the
process $X$ is supposed to be the Lévy driven Ornstein-Uhlenbeck process
defined as:

$dX_{t}=-X_{t}dt+2dZ_{t},$

where the driving noise $Z$ is the Lévy process with
$Z_{1}\sim\text{NIG}(1,0,1,0)$. The normal inverse Gaussian (NIG) random
variable is defined as the normal-mean variance mixture of the inverse
Gaussian random variable, and the probability density function of
$Z_{t}\sim\text{NIG}(\alpha,\beta,\delta t,\mu t)$ is given by

$x\mapsto\frac{\alpha\delta t\exp\\{\delta
t\sqrt{\alpha^{2}-\beta^{2}}+\beta(x-\mu t)\\}K_{1}(\alpha\psi(x;\delta t,\mu
t))}{\pi\psi(x;\delta t,\mu t)}$

where $\alpha^{2}:=\gamma^{2}+\beta^{2}$ and $\psi(x;\delta t,\mu
t):=\sqrt{(\delta t)^{2}+(x-\mu t)^{2}}$. More detailed theoretical properties
of the NIG-Lévy process are given for example in [2].

For the simulation by YUIMA , we formally introduce the following system:

(4.4) $d\Bigg{[}\begin{array}[]{c}X_{1,t}\\\
X_{2,t}\end{array}\Bigg{]}=\Bigg{[}\begin{array}[]{c}X_{2,t}\\\
-X_{2,t}\end{array}\Bigg{]}dt+\Bigg{[}\begin{array}[]{c}0\\\
2\end{array}\Bigg{]}dZ_{t},$

with the initial condition:

$X_{1,0}=0,\ X_{2,0}=0.$

The process $X_{1,t}$ corresponds to the regressor $\int_{0}^{t}X_{u}du$ in
the regression model (4.3). Due to the specification of the new YUIMA
function, we will additionally construct the “full” model:

(4.5) $Y_{t}=\mu_{1}X_{1,t}+\mu_{2}X_{2,t}+\sigma J_{t},\quad J_{1}\sim
t_{\nu},$

and simulate the trajectory of $Y$ with $\mu_{2}=0$. After that, we will
estimate the parameters $(\mu_{1},\sigma,\nu)$ by the original model (4.3).
From now on, we show the implementation of this model in YUIMA . The sampling
setting is unchanged from the previous example code:
$h_{n}=1/50,n=1/h_{n}=50,T_{n}=50,B_{n}=15$.

R> library(yuima)########### Inputs ############ Inputs for Fourier
InversionR> method_Fourier = "FFT"; up = 6; low = -6; N = 2^17; N_grid =
60000# Inputs for the sample grid timeR> Factor <\- 1R> Factor1 <\- 1R>
initial <\- 0; Final_Time <\- 50 * Factor; h <\- 0.02/Factor1

The role of the variable Factor1 is the same as in the previous section. By
using the constructor `setLaw_th`, we define yuima.th-class in a similar
manner. After that, we construct the “full” model (4.5) by the constructors
setModel and setLRM. A trajectory of $Y$ and its regressors can be simulated
by the YUIMA method simulate.

##################################################### Regressor: Integrated
NIG-Levy driven OU process
#####################################################R> mu1 <\- -3; mu2 <\- 0;
scale <\- 3; nu <\- 2.5 # Model Parameters# Model DefinitionR> lawILOU <\-
setLaw_th(method = method_Fourier, up = up, low = low, N = N,\+ N_grid =
N_grid) # yuima.thR> regrILOU <\- setModel(drift = c("X2", "-X2"), jump.coeff
= c("0", "2"),\+ solve.variable = c("X1","X2"), xinit = c(0,0), measure.type =
"code",\+ measure = list(df = "rNIG(z, 1, 0, 1, 0)")) # Regressors definition#
t-Regression modelR> ModILOU <\- setLRM(unit_Levy = lawILOU, yuima_regressors
= regrILOU)# SimulationR> sampILOU <\- setSampling(initial, Final_Time, n =
Final_Time/h)R> true.parILOU <\- unlist(list(mu1 = mu1, mu2 = mu2, sigma0 =
scale, nu = nu))R> set.seed(1)R> simILOU <\- simulate(ModILOU, true.parameter
= true.parILOU, sampling = sampILOU)

Next we extract the trajectories of $Y$ and $X_{1}$ from the yuima.LevyRM-
class: simILOU and construct the new model for the estimation of the
parameters. Figure 5 shows the simulated sample paths for the regressor and
the response process $Y$.

Figure 5. Simulated Trajectory of the Student Lévy Regression model defined in
(4.3).

# Data extractionR> Dataset<-get.zoo.data(simILOU)R> newData <-Dataset[-2]R>
newData <\- merge(newData$X1, newData$Y)R> colnames(newData) <\- c("X1","Y")R>
plot(newData)# Define the model for estimationR> regrILOU1 <\- setModel(drift
= c("0"), diffusion = matrix(c("0"),1,1),\+ solve.variable = c("X1"), xinit =
c(0)) # Regressors definition# t-Regression modelR> ModILOU1 <\-
setLRM(unit_Levy = lawILOU, yuima_regressors = regrILOU1)# EstimationR> lower1
<\- list(mu1 = -10, sigma0 = 0.1)R> upper1 <\- list(mu1 = 10, sigma0 =
10.01)R> startILOU1<\- list(mu1 = runif(1, -10, 10), sigma0 = runif(1, 0.01,
4))R> Bn <\- 15*FactorR> Data1 <\- setData(newData)R> estILOU1 <\-
estimation_LRM(start = startILOU1, model = ModILOU1, data = Data1,\+ upper =
upper1, lower = lower1, PT = Bn)R> summary(estILOU1)

Quasi-Maximum likelihood estimationCall:estimation_LRM(start = startILOU1,
model = ModILOU1, data = Data1, upper = upper1, lower = lower1, PT =
Bn)Coefficients: Estimate Std. Errormu1 -2.959621 0.02112734sigma0 2.778659
0.03208519nu 2.404843 0.13018410-2 log L: 69711.32 854.3838

#### Example 2

In this example, we consider the following regression model:

(4.6) $Y_{t}=\mu_{1}\int_{0}^{t}V_{1,u}du+\mu_{2}\int_{0}^{t}V_{2,u}du+\sigma
J_{t},\ \quad J_{1}\sim t_{\nu}$

where the process $X=(V_{1},V_{2})$ satisfy

(4.7) $\displaystyle
dV_{1,t}=\frac{1}{\epsilon}(V_{1,t}-V_{1,t}^{3}-V_{2,t}+s)dt,$ (4.8)
$\displaystyle dV_{2,t}=(\gamma
V_{1,t}-V_{2,t}+\alpha)dt+\sigma^{\prime}dw_{t},$

with

$(\epsilon,s,\gamma,\alpha,\sigma^{\prime})=(1/3,0,3/2,1/2,2),$

and standard Wiener process $w$. We set the true values

$(\mu_{1,0},\mu_{2,0},\sigma_{0},\nu_{0})=(8,-4,8,3).$

The process $V$ is the so-called stochastic FitzHugh-Nagumo process which is a
classical model for describing a neuron; $V_{1}$ expresses the membrane
potential of the neuron and $V_{2}$ represents a recovery variable. Its
theoretical properties such as hypoellipticity, and Feller and mixing
properties are well summarized in [17]. The paper also provides a
nonparametric estimator of the invariant density and spike rate. Similarly,
other integrated degenerate diffusion can be considered as regressors. Its
ergodicity for the regularity condtions is studied for example in [25]. For
the statistical inference for degenerate diffusion processes, we refer to [9]
and [12].

For the implementation of the regression model (4.6) on YUIMA, we formally
consider the following dynamics:

(4.9) $d\left[\begin{array}[]{c}V_{1,t}\\\ V_{2,t}\\\ V_{3,t}\\\
V_{4,t}\end{array}\right]=\left[\begin{array}[]{c}V_{3,t}\\\ V_{4,t}\\\
3(V_{3,t}-V_{3,t}^{3}-V_{4,t})\\\
\frac{3}{2}V_{3,t}-V_{4,t}+\frac{1}{2}\end{array}\right]dt+\left[\begin{array}[]{c}0\\\
0\\\ 0\\\ 2\end{array}\right]dw_{t},$

with the initial condition:

$V_{1,0}=0,\ V_{2,0}=0,\ V_{3,0}=0,\ V_{4,0}=0.$

The first and second elements correspond to the regressor in (4.6). As in the
previous example, we first simulate data by the “full” model defined as:

(4.10)
$Y_{t}=\mu_{1}V_{1,t}+\mu_{2}V_{2,t}+\mu_{3}V_{3,t}+\mu_{4}V_{4,t}+\sigma
J_{t},\quad J_{1}\sim t_{\nu},$

with $\mu_{3}=\mu_{4}=0$, and after that, we extract the simulated data and
estimate the parameters based on the original regression model (4.6). Below we
show how to implement on YUIMA. In this example, we set
$h_{n}=1/200,n=1/h_{n}=200,T_{n}=1000,B_{n}=300$. The values of them are
controlled by the variables Factor and Factor1 in the example code.

R> library(yuima)########### Inputs ############ Inputs for Fourier
InversionR> method_Fourier = "FFT"; up = 6; low = -6; N = 2^17; N_grid =
60000# Inputs for the sample grid timeR> Factor <\- 20R> Factor1 <\- 4R>
initial <\- 0; Final_Time <\- 50 * Factor; h <\-
0.02/Factor1#############################################################
Regressor: Integrated stochastic FitzHugh-Nagumo process
#############################################################R> mu1 <\- 8; mu2
<\- -4; mu3 <\- 0; mu4 <\- 0; scale <\- 8; nu <\- 3 # Model Parameters# Model
DefinitionR> lawFN <\- setLaw_th(method = method_Fourier, up = up, low =
low,\+ N = N, N_grid = N_grid) # yuima.thR> regrFN <\- setModel(drift =
c("V3", "V4", "3*(V3-V3^3-V4)", "1.5*V3-V4+0.5"),\+ diffusion = matrix(c("0",
"0", "0", "2"),4,1),\+ solve.variable = c("V1","V2", "V3", "V4"),\+ xinit =
c(0,0,0,0)) # Regressors definitionR> ModFN <\- setLRM(unit_Levy = lawFN,
yuima_regressors = regrFN) # t-regression model# SimulationR> sampFN <\-
setSampling(initial, Final_Time, n = Final_Time/h)R> true.parFN <\-
unlist(list(mu1 = mu1, mu2 = mu2, mu3 = mu3, mu4 = mu4,\+ sigma0 = scale, nu =
nu))R> set.seed(12)R> simFN <\- simulate(ModFN, true.parameter = true.parFN,
sampling = sampFN)

In this example, both Factor and Factor1 are larger than the previous example,
and hence it takes a higher computational load than the previous one. Figure 6
shows the simulated sample paths for the regressor and the Student Lévy
Regression model $Y$. Now we move on to the estimation phase.

Figure 6. Simulated Trajectory of the Student Lévy Regression model defined in
(4.6).

R> Dataset<-get.zoo.data(simFN)R> newData <-Dataset[-c(3,4)]R> newData <\-
merge(newData$X1,newData$X2, newData$Y)R> colnames(newData) <\- c("X1","X2",
"Y")R> plot(newData)# Define the model for estimationR> regrFN1 <\-
setModel(drift = c("0", "0"), diffusion = matrix(c("0", "0"),2,1),\+
solve.variable = c("X1", "X2"), xinit = c(0,0)) # Regressors definitionR>
ModFN1 <\- setLRM(unit_Levy = lawFN, yuima_regressors = regrFN1) #
t-Regression model.# EstimationR> lower1 <\- list(mu1 = -10, mu2 = -10, sigma0
= 0.1)R> upper1 <\- list(mu1 = 10, mu2 =10, sigma0 = 10.01)R> startFN1<\-
list(mu1 = runif(1, -10, 10), mu2 = runif(1, -10, 10),\+ sigma0 = runif(1,
0.01, 4))R> Bn <\- 15*FactorR> Data1 <\- setData(newData)R> estFN1 <\-
estimation_LRM(start = startFN1, model = ModFN1, data = Data1,\+ upper =
upper1, lower = lower1, PT = Bn)R> summary(estFN1)

Quasi-Maximum likelihood estimationCall:estimation_LRM(start = startFN1, model
= ModFN1, data = Data1, upper = upper1, lower = lower1, PT = Bn)Coefficients:
Estimate Std. Errormu1 8.041838 0.04895791mu2 -4.041518 0.04495255sigma0
7.712157 0.04452616nu 3.020972 0.11837331-2 log L: 403960.9 1291.516

###  4.3. Real data regressors

In this section, we show how to use the YUIMA package for the estimation of a
Student Lévy Regression model in a real dataset. We consider a model where the
daily price of the Standard and Poor 500 Index is explained by the VIX index
and two currency rates: the YEN-USD and the Euro Usd rates. The dataset was
provided by Yahoo.finance and ranges from December 14th 2014 to May 12th 2023.
We downloaded the data using the function getSymbol available in quantmod
library. To consider a small value for the step size $h$ we estimate our model
every month ($h=1/30$). We report below the code for storing the market data
in an object of yuima.data-class

R> library(yuima)R> library(quantmod)R> getSymbols("^SPX", from =
"2014-12-04", to = "2023-05-12")# RegressorsR> getSymbols("EURUSD=X", from =
"2014-12-04", to = "2023-05-12")R> getSymbols("VIX", from = "2014-12-04", to =
"2023-05-12")R> getSymbols("JPYUSD=X", from = "2014-12-04", to =
"2023-05-12")R> SP <\- zoo(x = SPX$SPX.Close, order.by =
index(SPX$SPX.Close))R> Vix <\- zoo(x = VIX$VIX.Close/1000, order.by =
index(VIX$VIX.Close))R> EURUSD <\- zoo(x = ‘EURUSD=X‘[,4], order.by =
index(‘EURUSD=X‘))R> JPYUSD <\- zoo(x = ‘JPYUSD=X‘[,4], order.by =
index(‘JPYUSD=X‘))R> Data <\- na.omit(na.approx(merge(Vix, EURUSD, JPYUSD,
SP)))R> colnames(Data) <\- c("VIX", "EURUSD", "JPYUSD", "SP")R> days <\-
as.numeric(index(Data))-as.numeric(index(Data))[1]# equally spaced grid time
dataR> Data <\- zoo(log(Data), order.by = days)R> Data_eq <\- na.approx(Data,
xout = days[1] : tail(days,1L))R> yData <\- setData(zoo(Data_eq, order.by =
index(Data_eq)/30)) # Data on monthly basis We decided to work with log-price
(see Data variable in the above code) to have quantities defined on the same
support of a Student-$t$ Lévy process, i.e.: the real line. Since the
estimation method in YUIMA requires that the data are observed on an equally
spaced grid time, we interpolate linearly to estimate possible missing data to
get a log price for each day. Figure 7 shows the trajectory for each financial
series.

Figure 7. Financial data observed on a monthly equally grid time ranging from
December 14th 2014 to May 12th 2023.

To estimate the model we perform the same steps discussed in the previous
examples and we report below the code for reproducing our result.

# Inputs for integration in the inversion formulaR> method_Fourier <\- "FFT";
up <\- 6; low <\- -6; N <\- 10000; N_grid <\- 60000# Model DefinitionR> law
<\- setLaw_th(method = method_Fourier, up = up, low = low, N = N,\+ N_grid =
N_grid)R> regr <\- setModel(drift = c("0", "0", "0"), diffusion =
matrix("0",3,1),\+ solve.variable = c("VIX", "EURUSD", "JPYUSD"), xinit =
c(0,0,0))Mod <\- setLRM(unit_Levy = law, yuima_regressors = regr, LevyRM =
"SP")# EstimationR> lower <\- list(mu1 = -100, mu2 = -200, mu3 = -100, sigma0
= 0.01)R> upper <\- list(mu1 = 100, mu2 = 100, mu3 = 100, sigma0 = 200.01)R>
start <\- list(mu1 = runif(1, -100, 100), mu2 = runif(1, -100, 100),\+ mu3 =
runif(1, -100, 100), sigma0 = runif(1, 0.01, 100))R> est <\-
estimation_LRM(start = start, model = Mod, data = yData, upper = upper,\+
lower = lower, PT = floor(tail(index(Data_eq)/30,1L)/2))R> summary(est)

Quasi-Maximum likelihood estimationCall:estimation_LRM(start = start, model =
Mod, data = yData, upper = upper, lower = lower, PT =
floor(tail(index(Data_eq)/30, 1L)/2))Coefficients: Estimate Std. Errormu1
0.000335328 0.004901559mu2 -0.062235188 0.033027232mu3 0.016291220
0.039382559sigma0 0.082031691 0.004859138nu 3.062728822 0.616466913-2 log L:
-1506.067 48.17592

##  5\. Conclusion

In this paper, we have presented classes and methods for a $t$-Lévy regression
model. In particular, the simulation and the estimation algorithm have been
introduced from a computational point of view. Moreover three different
algorithms have been developed for computing the density, the cumulative
distribution and the quantile functions and the random number generator of the
$t$-Lévy increments defined over a non-unitary interval. These latter methods
can be also used in any stochastic model available in YUIMA for the definition
of the underlying noise.

Based on our simulated data, for the estimation of the degrees of freedom the
Laguarre method seems to be more accurate. However, due to the restriction on
the number of roots used in the evaluation of the integral, we notice a bias
on the estimation of the scale parameter. Such bias seems to be observable
also in the other implemented methods when the same number of points is used
in the evaluation of the integral for the inversion of the characteristic
function. However the COS and the FFT methods, the numerical precision can be
improved. In particular, it is possible to ensure that the estimates fall in
the asymptotic confidence interval at 95% level even if the estimates of the
degrees of freedom in the simplest model (first example) seem to underestimate
the true value. As concerns computational time, the FFT methods with a
sufficiently large computational precision ($N\geq 5000$) and $h$ sufficiently
small seems to be an acceptable compromise.

## Acknowledgment

This work is supported by JST CREST Grant Number JPMJCR2115, Japan and by the
MUR-PRIN2022 Grant Number 20225PC98R, Italy CUP codes: H53D23002200006 and
G53D23001960006.

## References

  * [1] S. Asmussen and J. Rosiński.  Approximations of small jumps of lévy processes with a view towards simulation.  Journal of Applied Probability, 38(2):482–493, 2001.
  * [2] O. E. Barndorff-Nielsen.  Processes of normal inverse Gaussian type.  Finance Stoch., 2(1):41–68, 1998.
  * [3] C. Berg and C. Vignat.  On some results of cufaro petroni about student t-processes.  Journal of Physics A: Mathematical and Theoretical, 41(26):265004, 2008.
  * [4] A. Brouste, M. Fukasawa, H. Hino, S. Iacus, K. Kamatani, Y. Koike, H. Masuda, R. Nomura, T. Ogihara, Y. Shimuzu, M. Uchida, and N. Yoshida.  The yuima project: A computational framework for simulation and inference of stochastic differential equations.  Journal of Statistical Software, 57(4):1-51, 2014.
  * [5] L. V. Castro-Porras, R. Rojas-Martínez, C. A. Aguilar-Salinas, O. Y. Bello-Chavolla, C. Becerril-Gutierrez, and C. Escamilla-Nuñez.  Trends and age-period-cohort effects on hypertension mortality rates from 1998 to 2018 in mexico.  Scientific Reports, 11(1):17553, 2021.
  * [6] J. W. Cooley and J. W. Tukey.  An algorithm for the machine calculation of complex fourier series.  Mathematics of computation, 19(90):297–301, 1965.
  * [7] N. Cufaro Petroni.  Mixtures in nonstable Lévy processes.  J. Phys. A, 40(10):2227–2250, 2007.
  * [8] L. Devroye.  On the computer generation of random variables with a given characteristic function.  Comput. Math. Appl., 7(6):547–552, 1981.
  * [9] S. Ditlevsen and A. Samson.  Hypoelliptic diffusions: filtering and inference from complete and partial observations.  J. R. Stat. Soc. Ser. B. Stat. Methodol., 81(2):361–384, 2019.
  * [10] F. Fang and C. W. Oosterlee.  A novel pricing method for european options based on fourier-cosine series expansions.  SIAM Journal on Scientific Computing, 31(2):826–848, 2009.
  * [11] F. Fang and C. W. Oosterlee.  A fourier-based valuation method for bermudan and barrier options under heston’s model.  SIAM Journal on Financial Mathematics, 2(1):439–463, 2011.
  * [12] A. Gloter and N. Yoshida.  Adaptive estimation for degenerate diffusion processes.  Electron. J. Stat., 15(1):1424–1472, 2021.
  * [13] S. Haberman and A. Renshaw.  On age-period-cohort parametric mortality rate projections.  Insurance: Mathematics and Economics, 45(2):255–270, 2009.
  * [14] C. C. Heyde and N. N. Leonenko.  Student processes.  Adv. in Appl. Probab., 37(2):342–365, 2005.
  * [15] F. Hubalek.  On simulation from the marginal distribution of a student t and generalized hyperbolic lévy process.  2005\.
  * [16] S. M. Iacus and N. Yoshida.  Simulation and inference for stochastic processes with yuima.  A comprehensive R framework for SDEs and other stochastic processes. Use R, 2018.
  * [17] J. R. León and A. Samson.  Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate.  Ann. Appl. Probab., 28(4):2243–2274, 2018.
  * [18] A. Loregian, L. Mercuri, and E. Rroji.  Approximation of the variance gamma model with a finite mixture of normals.  Statistics & Probability Letters, 82(2):217–224, 2012.
  * [19] T. Massing.  Simulation of student–lévy processes using series representations.  Computational Statistics, 33(4):1649–1685, 2018.
  * [20] H. Masuda, L. Mercuri, and Y. Uehara.  Noise inference for ergodic Lévy driven SDE.  Electronic Journal of Statistics, 16(1):2432 – 2474, 2022.
  * [21] H. Masuda, L. Mercuri, and Y. Uehara.  Quasi-likelihood analysis for student-lévy regression.  arXiv preprint arXiv:2306.16790, 2023.
  * [22] L. Mercuri, A. Perchiazzo, and E. Rroji.  Finite mixture approximation of carma (p, q) models.  SIAM Journal on Financial Mathematics, 12(4):1416–1458, 2021.
  * [23] R. Singleton.  An algorithm for computing the mixed radix fast fourier transform.  IEEE Transactions on audio and electroacoustics, 17(2):93–103, 1969\.
  * [24] C. Sørensen.  Modeling seasonality in agricultural commodity futures.  Journal of Futures Markets: Futures, Options, and Other Derivative Products, 22(5):393–426, 2002.
  * [25] L. Wu.  Large and moderate deviations and exponential convergence for stochastic damping Hamiltonian systems.  Stochastic Process. Appl., 91(2):205–238, 2001.