# Predicting Stellar Mass Accretion: An Optimized Echo State Network Approach
in Time Series Modeling

Gianfranco Bino1 Shantanu Basu2,4 Ramit Dey2 Sayantan Auddy5 Lyle Muller3
Eduard I. Vorobyov6,7 1Department of Applied Mathematics, University of
Western Ontario, London, ON, N6A 5B7, Canada  2Department of Physics &
Astronomy, University of Western Ontario, London, ON, N6A 3K7, Canada
3Department of Mathematics, University of Western Ontario, London, ON, N6A
5B7, Canada  4Institute for Earth & Space Exploration, University of Western
Ontario, London, Ontario, N6A 5B7, Canada  5Jet Propulsion Laboratory,
California Institute of Technology, Pasadena, CA 91109, USA  6Department of
Astrophysics, The University of Vienna, A-1180 Vienna, Austria  7Research
Institute of Physics, Southern Federal University, Rostov-on-Don 344090,
Russia

###### Abstract

Modeling the dynamics of the formation and evolution of protostellar disks as
well as the history of stellar mass accretion typically involve the numerical
solution of complex systems of coupled differential equations. The resulting
mass accretion history of protostars is known to be highly episodic due to
recurrent instabilities and also exhibits short timescale flickering. By
leveraging the strong predictive abilities of neural networks, we extract some
of the critical temporal dynamics experienced during the mass accretion
including periods of instability. Particularly, we utilize a novel form of the
echo state neural network (ESN), which has been shown to deal efficiently with
data having inherent nonlinearity. We introduce the use of optimized-ESN (Opt-
ESN) to make model-independent time series forecasting of mass accretion rate
in the evolution of protostellar disks. We apply the network to multiple
hydrodynamic simulations with different initial conditions and exhibiting a
variety of temporal dynamics to demonstrate the predictability of the Opt-ESN
model. The model is trained on simulation data of $\sim 1-2$ Myr, and achieves
predictions with a low normalized mean square error ($\sim 10^{-5}$ to
$10^{-3}$) for forecasts ranging between 100 and 3800 yr. This result shows
the promise of the application of machine learning based models to time-domain
astronomy.

###### keywords:

Stellar accretion (1578) — Neural networks (1933) — Star formation (1569)

##  1 Introduction

We are entering a new era of rapid advance in time-domain astronomy that
promises to revolutionize our understanding of transient astrophysical
phenomena. These advances will occur through a variety of instruments that
span the electromagnetic and gravitational-wave spectrum. It is important to
push forward the development of analysis techniques that can in principle be
utilized to model all transient phenomena regardless of the signal source or
parameters.

In this paper, we focus on modeling luminosity variations in the evolution of
young stellar objects (YSOs). These are precursors of main-sequence stars that
form from the collapsing dense regions of interstellar molecular clouds. Stars
in their early stages of evolution accumulate materials via mass accretion
from the surrounding accretion disk. Mass accretion from the disk to the
central object in the early evolutionary stage is likely driven by
gravitational torques arising from nonaxisymmetric spiral waves (Vorobyov &
Basu, 2007, 2009). Other drivers such as disk winds (Bai, 2013; Suzuki et al.,
2016), hydrodynamic and/or magnetohydrodynamic turbulence (Balbus, 2003), etc.
also lead to redistribution of angular momentum, enabling accretion of disk
material onto the central star. However, continuing mass infall onto the disk
from larger distances often leads to sustained or recurring gravitational
instability (GI) (Vorobyov & Basu, 2005, 2006) in the disk. The GI further
triggers the formation of gas clumps that move inward, resulting in bursts of
mass accretion onto the central object.

Such episodic accretion onto YSOs is frequently observed and is well known in
the form of FU Orionis objects (FUors) and EX Lupi objects (EXors). Recently,
episodic accretion bursts have also been detected in young massive protostars
(Caratti o Garatti et al., 2017). The FUors show a rapid rise of luminosity
from a few $L_{\odot}$ to $100-300$ $L_{\odot}$ (Audard et al., 2014). This
typically corresponds to an increase in mass accretion rate from $10^{-7}$
$M_{\odot}$ yr-1 to a few times $10^{-4}$ $M_{\odot}$ yr-1. The subsequent
decline of luminosity after the initial burst occurs over a timescale of many
decades. Due to the long timescale of decline, no FUor has ever been observed
to have more than one burst. In contrast, the EXors exhibit smaller luminous
amplitudes (up to a few tens $L_{\odot}$) in repetitive outbursts with
durations of several months. It is still uncertain whether these two phenomena
are related and if FUors correspond to an early stage of evolution with EXors
representing smaller amplitude bursts at a later stage (e.g., Contreras Peña
et al., 2019).

Furthermore, the mass accretion rates of YSOs are highly episodic due to
recurrent instabilities. They exhibit short timescale flickering due to
inherent nonlinearity and inhomogeneity in the disk structure (Elbakyan et
al., 2016). This makes forecasting burst events particularly challenging.
Developing analysis techniques using present-day simulation data is key to
advancing the study of such observations even if the underlying dynamics are
not well known.

The last decade has seen phenomenal growth in adaptations of various machine
learning (ML) techniques in analyzing astronomical data (Auddy et al., 2021,
2022) and making predictions in time-domain astronomy (Bloom & Richards, 2012;
Rocha-Solache et al., 2022). Neural network (NN) based models are particularly
powerful as they are not tied to a specific set of physical equations and
assumptions. They can be trained on data (both from simulation and
observation) to capture the nonlinear physics of the system and to make
predictions (Auddy & Lin, 2020). The objective of this paper is to demonstrate
that NN-based models can be used to forecast the evolution of transient
phenomena in real-time.

We introduce the use of an echo state neural network (ESN) (Lukoševičius,
2012; Kim & King, 2020) to make robust predictions of stellar mass accretion
of evolving YSOs. The model is trained on time-series data obtained from
hydrodynamical simulations (see for example Vorobyov & Basu, 2010) which
capture the evolution of such complex nonlinear star-disk systems. A series of
simulations (Vorobyov & Basu, 2005, 2006, 2010, 2015; Meyer et al., 2017) have
demonstrated the prevalence of such episodic accretion driven by mass infall
onto a nascent protostellar disk. In order to deal with the nonlinearity we
use a novel approach of dividing the (simulation) data into a slowly-varying
(“deterministic”) component and a more rapidly-varying (“fluctuating” or
“chaotic”) component. We train the ESN-based model on each component of the
data to make the subsequent prediction of the burst events. This ESN-based
framework lays the foundation for analyzing such transient phenomena from
upcoming surveys, like wide-field optical wavelength mapping with frequent
time sampling by the Zwicky Transient Facility (ZTF) and Vera C. Rubin
Observatory (VRO).

This paper is organized as follows. In Section 2 we discuss the hydrodynamic
simulations that capture the mass accretion in disk evolution. Section 3 gives
an overview of the ESN architecture. In Section 4 we introduce the Opt-ESN
model and outline the data preparation procedure. Results are presented in
Section 5. A further discussion is in Section 6 and conclusions are in Section
7. A reader who is focused on the astrophysical consequences of this work can
choose to read Section 2 and then move ahead to Section 5.

##  2 Hydrodynamic Simulations

Numerical simulations of disk evolution can be done using a set of
hydrodynamic equations that are vertically integrated along the direction of
the rotation axis and follow the nonaxisymmetric evolution of physical
variables in polar ($r,\phi$) coordinates. This is viable in the expected
scenario where the disk vertical scale height is significantly less than its
radial extent. A series of papers have employed the thin-disk approximation to
model the long term evolution of protostellar disks over several Myr
timescales (e.g., Vorobyov & Basu, 2006, 2010, 2015; Vorobyov et al., 2017,
2020). It is still challenging to model the full temporal range of disk
evolution using three-dimensional simulations, and state-of-the-art models
that resolve the central protostar can advance as far as $\sim 10^{3}$ yr past
protostar formation (see Machida & Basu, 2019)

We train the ESN on the long-term ($\sim 10^{6}$ yr) disk simulations
presented by Vorobyov et al. (2017). The simulations calculate the self-
consistent disk formation and evolution. This is done by starting from the
hydrodynamic collapse of a prestellar cloud core and continuing into the
protostellar phase with a central protostar and surrounding disk. The basic
equations and numerical finite difference numerical methods are described in
Vorobyov & Basu (2010) and Vorobyov et al. (2017). A numerical solution is
found to the partial differential equations describing the time and space
evolution of the mass surface density, the planar momentum components, and the
internal energy per unit area. Additional equations are employed to calculate
self-gravity, viscosity, and heating and cooling rates due to multiple
processes. A central sink cell of radius 5 au is adopted at the coordinate
origin in order to avoid very small time steps imposed by the Courant-
Friedrichs-Lewy condition, so that the long-term evolution of the remaining
region (radius $\sim 10^{4}$ au) can be calculated.

The solution of the disk evolution after protostar formation consists of a
highly episodic accretion process. While some features of the episodes can be
understood in a deterministic manner using the criterion for gravitational
instability (Das & Basu, 2022), the nonaxisymmetry and nonlinearity of the
problem lead to a time evolution of accretion rate that has stochastic and
chaotic features. Each simulation is quite costly in terms of run time (up to
several months on a single computer node with 48 physical cores), and a
typical parameter survey consists of $\sim 10$ models. What we explore here is
the possibility of taking a set of simulation models with different initial
conditions as input and training a neural network on a portion of the time
evolution in order to extract some intrinsic and underlying dynamics of the
system. We can then see how far in time the neural network can forecast the
solution into a regime where it was not trained.

###  2.1 Hydrodynamic Simulation Outputs

We utilize the simulation outputs in 6 of the 35 models presented in Vorobyov
et al. (2017). These six models differ in their initial conditions, which are
described here. The initial axisymmetric radial profiles for the gas surface
density $\Sigma$ and angular velocity $\Omega$ for the initial prestellar
collapsing core are

$\displaystyle\Sigma$
$\displaystyle=\frac{r_{0}\Sigma_{0}}{\sqrt{r^{2}+r_{0}^{2}}}\,,$ (1)
$\displaystyle\Omega$
$\displaystyle=2\,\Omega_{0}\left(\frac{r_{0}}{r}\right)^{2}\left[\sqrt{1+\left(\frac{r}{r_{0}}\right)^{2}}-1\right]\,,$
(2)

where $\Sigma_{0}$ and $\Omega_{0}$ are the surface mass density and angular
velocity at the center of the core, respectively. These are power-law profiles
with asymptotic dependence $\propto r^{-1}$ and have a central plateau radius
$r_{0}$ that is the length scale over which thermal pressure can smoothen the
density profile (for details, see Vorobyov et al., 2017).

To generate a gravitationally unstable core, each model is characterized by
the ratio $r_{\text{out}}/r_{0}=6$, where $r_{\text{out}}$ is the core’s outer
radius. The cloud core mass $M_{\text{cl}}$ is found using the initial radial
profile for the gas surface density $\Sigma$. The quantity $\Omega_{0}$ is
selected such that the models have an intial ratio of rotational to
gravitational energy $\beta_{0}$ in the range of $\approx 10^{-4}$ to $0.07$.
We summarize the initial model conditions in Table 1.

The simulation outputs of the mass accretion rate to the central sink
$\dot{M}(t)$ are shown in Figure 1 for each model. Comparison with the Table 1
parameter values shows that there is a general increase of the variability
amplitude as $M_{\rm core}$ and/or $\beta_{0}$ increase. Increasing mass or
angular momentum leads to more massive protostellar disks and greater activity
of GI induced bursts.

Figure 1: Mass accretion rate evolution from hydrodynamic simulations for each
of the six models labeled 26 through 32 (with the exception of model 28) in
Vorobyov et al. (2017). The mass accretion rate is shown in units of
$M_{\odot}$ yr-1.

##  3 Echo state neural networks

Neural networks (NN) have demonstrated the capability to approximate
continuous functions and are often referred to as universal approximators
(Schäfer & Zimmermann, 2006). In the context of time series analysis, this
gives them the ability to estimate the underlying dynamical processes
governing the system. This is done by using the NN as a mapping function
between the inputs and targeted outputs, allowing them to extract complex
temporal relationships within the time-series data. The NN architecture is
based on a collection of interconnected nodes, or “neurons”, as shown
schematically in Figure 2. The nodes are often arranged in layers from input
to output, as shown in Figure 3. These architectures can be further extended
to include recurrent units that maintain a network’s hidden state during model
training. These hidden states allow the network to recognize temporal
sequences in data, which can make a recurrent neural network (RNN)
(Salehinejad et al., 2017) particularly useful in the context of time series
analysis.

Figure 2: Schematic diagram of a single layer perceptron. The input data and
bias are given as $\textbf{x}=[b,x_{1},x_{2},\ldots,x_{n}]$ while the weights
are given as $\textbf{W}=[w_{b},w_{1},w_{2},\ldots,w_{n}]$. $f_{\text{act}}$
is the activation function through which $\sum_{i}W_{i}x_{i}$ is passed to
return an output $y=f_{act}\left(\sum_{i}W_{i}x_{i}\right)$.  Figure 3:
Schematic diagram of a multilayer perceptron neural network with several
hidden layers. The neurons in the hidden layers take a linear combination of
the inputs from the previous layer and passes it through an activation
function to generate an output, which is further passed to the set of neurons
in the next layer. Figure 4: General architecture of an echo state neural
network, shown with sparsely connected random reservoir units. Here,
$\textbf{W}_{\text{input}}$ and $\textbf{W}_{r}$ are randomly generated sparse
matrices.  Figure 5: Input to output sequence for both the validation (left)
and testing (right) sets. The validation set is considered “in-time” because
each output value $y_{i}$ is calculated using an input $x_{i}$ taken directly
from the given data set. The testing set however, is considered “out-of-time”
because the model output is recursively served as an input to the following
time step. That is, given the data set’s terminal point $x_{n}$, the
prediction of $x(t)$ at the following point in time is given by the network
output $y_{n}$. As we step beyond the horizon of the given data set, $y_{n}$
will be used as the input to calculate the prediction $y_{n+1}$ and so on.

An ESN can be considered to be a sparsely connected RNN where the hidden
layers along with the weights act as a “reservoir”. This reservoir functions
as a nonlinear temporal kernel, embedding the dynamics of the input data onto
a higher dimensional computation space. For an ESN architecture only the
reservoir-to-output weights are trainable while the input-to-hidden and
hidden-to-hidden weights are chosen randomly and kept fixed during the
training process. The sparsity of the ESN architecture and the fact that the
hidden layer weights are not updated during the training process,
automatically addresses the problem of vanishing gradients, as seen typically
for a more conventional RNN based model (Jaeger, 2007; Lukoševičius, 2012).
Since only the output weights are trainable, which is a simple linear
regression task compared to the slow convergence of tuning the parameters of
other networks, ESNs are much faster to train compared to other RNNs. For
chaotic time-series prediction, ESNs have shown exceptionally good performance
as these networks can capture the nonlinear dynamics of the system
efficiently.

In order to effectively model the chaotic dynamics that govern the
hydrodynamic simulations, we implement the use of an ESN based model. An ESN
architecture having a input, reservoir and output layer with its corresponding
weights is demonstrated in Figure 4. For an input time series $\textbf{x}(t)$,
we begin by defining the following input and target time sequences:

$\displaystyle\textbf{x}_{1}(t)=[x_{1}(t),x_{2}(t),\cdots,x_{n-1}(t)]\quad\;\text{Input
Sequence}$
$\displaystyle\textbf{x}_{2}(t)=[x_{2}(t),x_{3}(t),\cdots,x_{n}(t)]\quad\quad\;\text{Target
Sequence}$

In order to utilize the neural network $\mathcal{N}$ as a predictive time
series model, we form the mapping
$\mathcal{N}:\textbf{x}_{1}(t)\mapsto\textbf{x}_{2}(t)$ to extract the
relationship between the quantities $x_{i}(t)$ and $x_{i+1}(t)$. To do so, we
train the ESN using the following steps:

  * •

Randomly generate the the input weight matrix $\textbf{W}_{\text{input}}$ and
the reservoir weight matrix $\textbf{W}_{r}$.

  * •

For each quantity $x_{i}(t)$ in $\textbf{x}_{1}(t)$, construct an $N_{r}\times
1$ reservoir state vector $\textbf{v}_{i}$, initialized to
$\textbf{v}_{1}=\textbf{0}$. Let

$\displaystyle\textbf{v}_{i+1}=(1-\alpha)\textbf{v}_{i}+\alpha
f_{\mathrm{act}}(\mathcal{W}_{i})\,,$ (3)
$\displaystyle\text{with}\quad\mathcal{W}_{i}=\textbf{W}_{\mathrm{input}}x_{i}(t)+\textbf{W}_{r}\textbf{v}_{i}+\textbf{W}_{b}\,,$

where $N_{r}$ is the reservoir size and $0<\alpha<1$ is the leaking rate.
Equation (3) includes a randomly generated bias term $\textbf{W}_{b}$ and
adopts the activation function $f_{\mathrm{act}}(x)=\tanh(x)$.

  * •

Define a washout quantity $\omega<n$ as an initially discarded transient and
for every $i>\omega$, construct the internal state

$\textbf{X}=\begin{bmatrix}1&1&\cdots&1\\\
x_{i}(t)&x_{i+1}(t)&\cdots&x_{n-1}(t)\\\
\textbf{v}_{i}&\textbf{v}_{i+1}&\cdots&\textbf{v}_{n-1}\\\ \end{bmatrix}.$ (4)

  * •

Finally, compute the output matrix using the Moore-Penrose inverse on the set
$\overline{\textbf{x}}_{2}(t)=[x_{\omega+2}(t),x_{\omega+3}(t),\cdots,x_{n}(t)]$,
yielding

$\textbf{W}_{\mathrm{output}}=\overline{\textbf{x}}_{2}^{\ast}\left(\textbf{X}^{\ast}\textbf{X}\right)^{-1}\textbf{X}^{\ast}.$
(5)

Note that if the matrix $\left(\textbf{X}^{\ast}\textbf{X}\right)$ is near-
singular, it is recommended to regularize the regression by adding a constant
$\lambda$ along the diagonals (known as Tikhonov regularization).

Once $\textbf{W}_{\mathrm{output}}$ has been calculated, the output can be
computed as

$\textbf{y}(t)=\textbf{W}_{\mathrm{output}}\textbf{X}\,.$ (6)

This effectively represents an estimate to the mapping of one point in time to
the next. Therefore, in order to predict the $(i+1)$th time step (i.e.,
estimate the quantity $x_{i+1}$), we construct the internal state v with
$x_{i}$ and use $\textbf{W}_{\mathrm{output}}$ to compute the output using
Equation (6).

The hyperparameters used for defining the reservoir and characterzing the
network are described below:

  * •

The reservoir size $N_{r}$
Determines the number of units in the reservoir (or in turn the size of the
reservoir).

  * •

Spectral Radius $\rho$
This is a global parameter that determines the maximal eigenvalue of the
$\textbf{W}_{r}$ matrix. In other words it scales the reservoir connection
matrix and controls the width of the distribution of the nonzero elements
present in $\textbf{W}_{r}$. In most cases, $\rho(\textbf{W})<1$ maintains the
echo state property.

  * •

Input scaling $\varrho$
This parameter determines the scaling of the input weight matrix. It also
controls the amount of nonlinearity in the dynamics of the reservoir.

  * •

Connectivity $c_{r}$
Controls the degree of sparsity in the reservoir weight matrix.

  * •

Leaking Rate $\alpha$
Controls the speed of the reservoir dynamics in reaction to the input.

##  4 Optimized Echo State Neural Networks

###  4.1 Network Architecture

Liu et al. (2018) introduced a parallel series approach where one stacks a
series of reservoirs by generating $L$ independent input and reservoir
matrices. The time series is trained and validated through each reservoir to
form $L$ output matrices, Finally, the model’s output $\hat{\textbf{y}}(t)$ is
taken to be the mean of all $L$ realizations so that

$\hat{\textbf{y}}(t)=\frac{1}{L}\sum_{j}\textbf{y}^{(j)}(t)\,.$ (7)

The optimized-ESN (Opt-ESN) extends Equation (7) to being a weighted sum
rather than the standard mean. That is, the output realizations from each
reservoir are weighted by a set of optimal coefficients. As such, the final
output of the Opt-ESN is given by

$\hat{\textbf{y}}(t;\hat{\boldsymbol{\beta}})=\sum_{j}\hat{\beta}_{j}\textbf{y}^{(j)}(t)\,,$
(8)

where the coefficients $\hat{\boldsymbol{\beta}}$ are found by minimizing the
squared residuals over the input’s validation segment
$\textbf{x}_{\text{val}}(t)$. This is done by solving the linear optimization
problem to find the minimum value of the loss function

$\begin{array}[]{rl}\mathcal{L}=\left\lVert\textbf{x}_{\text{val}}(t)-\hat{\textbf{y}}(t;\hat{\boldsymbol{\beta}})\right\rVert^{2}\par\end{array}.$
(9)

Here, the validation segment is defined as the in-time portion of the data set
used to validate the model’s output (see Section 4.3.1 for more details).

###  4.2 Data Preparation

As the input to the Opt-ESN, we used a portion of the simulated data for each
model, corresponding to times of vigorous episodic accretion. We divide the
data into segments of various lengths having different time steps and use this
for training, validation and testing. Furthermore, as an aid to assessing the
quality of our forecasts, we characterize the time scale in terms of the
Lyapunov exponent. We define the dimensionless time length $\Lambda\cdot
N_{t}$ as the Lyapunov time, where $N_{t}$ is the observation number. Here,
the quantity $\Lambda$ is taken to be the maximum Lyapunov exponent which
characterizes the rate of separation between close trajectories in phase space
and effectively quantifies the degree of chaos present. That is, if two
trajectories are initially separated by some infinitesimal amount
$\Delta_{0}$, then the rate of divergence as a function of time $t$ is
approximately

$\displaystyle|\Delta(t)|\approx|\Delta_{0}|\cdot\text{exp}(\Lambda t).$ (10)

It becomes clear that for $\Lambda>0$, the separation $|\Delta(t)|$ grows
exponentially with time (Vulpiani et al., 2009). We estimate $\Lambda$ using
the algorithm in Eckmann et al. (1986) where our outputs are given in the
fourth column of Table 2. Our calculations demonstrate that each value of
$\Lambda$ is greater than zero, indicating that each of the simulation models
can be considered to be a chaotic system. We preprocess the simulation data
$\dot{M}(t)$ by assuming that it is separable in the form

$\dot{M}(t)=\dot{M}_{d}(t)+\dot{\mathcal{M}}(t)\,,$ (11)

where $\dot{M}_{d}(t)$ and $\dot{\mathcal{M}}(t)$ represent the data’s
deterministic and fluctuating components, respectively. In order to extract
the fluctuating component, we pass $\dot{M}(t)$ through a high pass
filter111We utilized Matlab’s highpass function with passband frequency fpass
= 2000 and sampling rate f = 50000. The units are the inverse observation
number.. The deterministic component can then be extracted by subtracting
$\dot{\mathcal{M}}(t)$ from $\dot{M}(t)$. Furthermore, we normalize each
component with respect to the standard deviation $\sigma$ of $\dot{M}(t)$. In
order to make predictions, we feed $\dot{M}_{d}(t)$ and $\dot{\mathcal{M}}(t)$
into the Opt-ESN separately, run the forecasts and sum the outputs to get the
final prediction on $\dot{M}(t)$. We found that decomposing and processing the
simulation data in this fashion gives better performing output than inputting
$\dot{M}(t)$ directly. The proposed network architecture is given in Figure 6.
The network has two important layers contributing to the output. At the
stacked reservoir layer, the output weight matrices
$\textbf{W}_{\text{output}}$ are determined using regression techniques over
the training segment. The weight matrices are then used to output a series of
unique paths over the validation segment where the coefficients
$\hat{\beta}_{j}$ are found by solving a linear optimization problem.

Figure 6: Schematic figure representing the proposed architecture of the Opt-
ESN. The input is separated into fluctuating and deterministic components and
fed through an Opt-ESN. Each of the individual forecasts are summed to produce
the estimated forecast to $\dot{M}(t)$. Here, the neural network is shown to
contain $L$ reservoirs stacked into a single layer used in formulating the
linear programming problem. Each reservoir is taken to be sparse and randomly
connected.

###  4.3 Hyperparameter Selection

Selecting the optimal set of hyperparameters can be difficult given the size
of the hyperparameter space. As such, for each component of the individual
simulation data set (fluctuating and deterministic), we perform a hybrid,
discrete-stochastic search over the entire hyperparameter space. That is, we
prespecify the search space for each hyperparameter, construct a nested loop
searching over all possible combinations, build a model and compute the mean
square error (MSE) on the validation segment. In addition to this, within each
iteration of the loop, we generate a random number for each hyperparameter
that is within the search space, construct and validate a second model, and
compute the respective MSE. These two models are then compared and the one
with the lowest MSE is kept. This process continues until the entire
hyperparameter space has been searched and the model with the lowest MSE on
the validation segment is selected.

Figure 7: The Opt-ESN predictions of the simulation data. The solid black line
shows a fraction of the simulation data used for training/validation, the
dotted black is the out-of-time testing segment of the data and the red line
is the network’s prediction.  Figure 8: Performance assessments on the
predictions made for the simulation data. The upper portion of the figure
demonstrates the rolling NMSE as a function of Lyapunov time for $\dot{M}(t)$,
$\dot{\mathcal{M}}(t)$ and $\dot{M}_{d}(t)$. The lower portion of the figure
gives the residual autocorrelation function on the validation segment.

####  4.3.1 Training and Data Segmentation

In order to train the Opt-ESN model, the input data must be split into three
segments:

  * •

Training Segment ($t_{0}<t\leq t_{\text{train}}$)
Used to train the model and compute the output matrix
$\textbf{W}_{\mathrm{output}}$.

  * •

Validation Segment ($t_{\text{train}}<t\leq t_{\text{val}}$)
Used as an in-time assessment on how well $\textbf{W}_{\mathrm{output}}$ maps
a single input to an output.

  * •

Testing Segment ($t_{\text{val}}<t\leq T$)
Used as a blind test for the model to assess how well the out-of-time
predictions perform.

The difference between the validation and testing segment is that the model is
dynamic within the validation segment. That is, each output is calculated
using an input that is directly from the validation segment of the data. In
the testing segment however, the model is using it’s own output as an input
for the next time step. This is demonstrated in Figure 5 where the validation
use is given on the left side diagram and the testing use is given on the
right side diagram.

Table 1: Initial conditions for each of the seven simulation models. Model | $M_{\text{core}}$ | $\beta_{0}$ | $r_{0}$ | $M_{*,\text{fin}}$
---|---|---|---|---
Model 26 | $1.245M_{\odot}$ | $1.27\%$ |  $2777$ au | $0.753M_{\odot}$
Model 27 | $1.076M_{\odot}$ | $0.56\%$ |  $2400$ au | $0.801M_{\odot}$
Model 29 | $0.999M_{\odot}$ | $0.28\%$ |  $2229$ au | $0.818M_{\odot}$
Model 30 | $1.537M_{\odot}$ | $1.27\%$ |  $3429$ au | $0.887M_{\odot}$
Model 31 | $1.306M_{\odot}$ | $0.28\%$ |  $2915$ au | $1.031M_{\odot}$
Model 32 | $1.383M_{\odot}$ | $0.56\%$ |  $3086$ au | $1.070M_{\odot}$

##  5 Results

###  5.1 Opt-ESN Outputs

In order to train the model, the simulation data for each model is decomposed
into deterministic and fluctuating components and standardized as discussed in
Section 4.2. Forecasting either component of the simulation data requires the
time series to be stationary. Nonstationary dynamics risk consequences such as
spurious correlations and heavily biased mean and variance estimates (see
Appendix for more details). We assess the stationarity of each component using
the augmented Dickey–Fuller (ADF) test at 5% significance (Patterson, 2011).
Ultimately, we find that in each model, $\dot{\mathcal{M}}(t)$ follows a
stationary process while $\dot{M}_{d}(t)$ does not. As such, we applied first
order differencing on the deterministic component to enforce stationarity.
Network hyperparameters are selected according to the methodology introduced
in Section 4.3. We utilize a stack of $L=10$ reservoirs for size $N_{r}=250$
neurons and solve the optimization problem for Equation (9) to generate the
final output. In Figure 7, the Opt-ESN forecast is given for all six
simulation models. The outputs are found by summing the individual forecasts
of $\dot{M}_{d}(t)$ and $\dot{\mathcal{M}}(t)$ (the individual forecasts are
shown in the Appendix).

We assess the Opt-ESN performance using the dimensionless normalized mean
square error

$\displaystyle\text{NMSE}=\frac{1}{n}\sum_{i=1}^{n}\frac{(y_{i}-\hat{y_{i}})^{2}}{\text{max}(\hat{y})-\text{min}(\hat{y})}\,,$
(12)

where $y_{i}$ and $\hat{y_{i}}$ are the observed and predicted values,
respectively, for $n$ data points. Here, we take $y=\dot{M}/\sigma$, which is
a dimensionless quantity. Lower NMSE values indicate stronger performance and
NMSE $=0$ equates to perfect accuracy. We consider any model having NMSE
$<10^{-2}$ to indicate good performance. The performance assessments are
summarized in Table 2 in the Appendix. Additionally, in the upper half of
Figure 8, the rolling NMSE on each model component is given as a function of
Lyapunov time. The Lyapunov time is the characteristic timescale on which a
system is chaotic and typically limits the predictability of the said system.
Here, the Lyapunov time is calculated using the maximum Lyapunov exponent,
$\Lambda=\Lambda_{\textrm{max}}$. If the model is well specified, the
residuals between the observed and predicted values are expected to be
approximately normally distributed and thus largely attributed to white noise.

Over the validation set, our goodness-of-fit assessment on the residuals
utilize the autocorrelation function. In the lower portion of Figure 8, the
autocorrelations are given for lagged residuals adopting a 99% and 95%
confidence interval represented by the solid boundary lines and dashed
boundary lines, respectively. Given that at least 90% of lagged residuals lie
within the confidence interval, we assume the process to be approximately
white noise. Furthermore, the residuals over the training phase are used to
assess the goodness-of-fit of the Opt-ESN. We achieve this using the one-
sample $t$-test for the null hypothesis that the residuals are sampled from a
normal distribution with zero mean. The test statistic is

$t^{*}=\frac{\bar{x}-\mu}{\sigma_{s}/\sqrt{n}}\,,$

where $\bar{x}$ is the sample mean, $\mu$ is the hypothesized mean,
$\sigma_{s}$ is the sample standard deviation and $n$ is the sample size. The
statistic follows a $t$-distribution with $n-1$ degrees of freedom. Our
results show that for each model, we cannot reject the null hypothesis at the
5% significance, indicating that the residuals are approximately normal with
zero mean. Figure 9 shows the residual distributions for each model.

Figure 9: Distribution of residuals for each simulation model after network
training. Goodness-of-fit assessed using the one-sample $t$-test for the null
hypothesis that the residuals are sampled from a normal distribution with mean
equal to zero.

####  5.1.1 Episodic Bursts and Stability Assessment

The simulation data across each model has demonstrated strong episodic
behavior where the dynamics exhibit multiple high magnitude bursts that are
driven largely by mass infall. The Opt-ESN has shown strong forecasting
capabilities on segments of the data where episodic bursts have not occurred.
An interesting hypothesis to test is whether the neural network can adequately
resolve the occurrence of a burst over a short time interval. To do so, we
follow the same protocol for data preparation and hyperparameter selection. We
isolate the occurrence of the first large burst in Model 27 and produce
forecasts over 50 time steps ($\approx 0.31$ kyr). To reflect some of the
challenges in modeling observational data, we train the Opt-ESN under various
conditions:

  1. 1.

Training with 10000 data points.

  2. 2.

Training with 5000, 2500 and 1000 data points.

  3. 3.

Training with 10000 data points with 1% and 5% noise added.

These conditions additionally serve to assess the stability in the model’s
predictive power. That is, under nonideal training conditions, we aim to
demonstrate that the network can still forecast meaningful outputs. Our
training conditions vary in severity including scenarios with added noise and
lack of data availability. We enumerate the training conditions below.
Training condition (1) represents an ideal scenario with no noise and the
entire dataset and has an NMSE of $8.29\times 10^{-3}$. Condition (2) varies
the training length, and condition (3) includes a moderate and high degree of
noise in the data. The Opt-ESN performance and summary for each condition is
given in Table 3 and demonstrated in Figure 11. It is evident that the Opt-ESN
demonstrates an ability to resolve the presence of an episodic burst with
noisy training data. With increased noise the magnitude suffers and the NMSE
increases relative to condition (1) by $\approx 37\%\,\rm{and}\,88\%$ when we
included 1% and 5% noise, respectively. Condition (2) gives insight into the
importance of data availability with respect to model performance. In this
scenario, compared to condition (1), the NMSE increases by $\approx 10\%$ and
$40\%$ when the training length is reduced to 5000 and 2500 data points,
respectively.

Figure 10: Episodic burst forecasts under noisy training conditions.
Conditions range in severity from no noise (left panel) to 5% noise addition
(right panel). Here, the degree of noise is taken as a fraction (1 or 5
percent) of the maximum value of mass accretion rate. The performance results
are given from left to right as $8.29\times 10^{-3}$, $1.84\times 10^{-2}$,
and $1.47\times 10^{-1}$, respectively. This demonstrates that the model
trained under condition (1) is performing the best.

Figure 11: Comparison of burst forecasts across all training conditions
(excluding 1000 training length for visualization purposes). Condition (1)
gives the lowest NMSE ($8.29\times 10^{-3}$), and is highlighted with
increased line width.

##  6 Discussion

As we enter a new era in time domain astronomy, leveraging robust predictive
models to make meaningful data inference is increasingly valuable. Neural
networks became an attractive class of algorithms that can be used to model
nonlinear data due to the universal approximation theorem (Cybenko, 1989;
Schäfer & Zimmermann, 2006). In particular, the reservoir computing framework
that governs the ESN is ideal in modeling time series that possess chaotic
(nonlinear) temporal structure. We used the Opt-ESN model to predict/forecast
the protostellar mass accretion rate by training it on simulated
hydrodynamical time series data. To avoid nonstationary dynamics, first-order
differencing was used on the deterministic components of the data. Our
methodology has generated robust and accurate outputs over several sets of
chaotic data (i.e., having a high degree of nonlinear temporal dynamics).

###  6.1 Short-term Burst Predictions and Data Availability

In Section 5.1.1, we utilized the Opt-ESN to predict an episodic burst under
various conditions. Our findings demonstrate that the main driver in model
degradation is a lack of data. Even with a relatively high degree of noise,
the network was able to resolve the presence of a burst within a 50 time step
forecast. The Opt-ESN forecasts are aimed to provide a introductory framework
as to how ML based models can be used to make inference on burst occurrences
and forecast mass accretion. This can possibly help in establishing any
potential relation between FUor and EXor phenomena. In practice, however,
having thousands of years of observational data on a single object is not
feasible. We see in Figure 11 that the model maintains relatively strong
predictive performance in scenarios with training lengths as low as 2500 data
points ($0.3092$ kyr). However, below this number, we see a significant
decrease in model performance. Therefore, addressing data availability becomes
a critical component in developing a robust neural network framework. In these
scenarios, we propose two approaches:

  1. 1.

Simulating synthetic data.
As demonstrated in our approach, we can leverage simulation data to train
neural networks in conjunction with observational data. Here, the simulated
data would be integrated with the observational data as part of the
training/validating segments of the model. The simulated data however should
reflect the physical properties of the observed object and should include
several occurrences of episodic bursts to assure proper coverage in the
possible dynamics.

  2. 2.

Aggregating data from similar objects.
An alternative to generating synthetic data can be to aggregate data across
several observed objects. That is, the neural network is trained across
multiple similar objects. An advantage in this approach is that no portion in
the training data is synthetic and as such, the forecasts will reflect the
true dynamics of the observations. This approach can alleviate some of the
data availability issues in practice, however it is crucial that each of the
observed objects reflect the physical properties of the entire system being
modeled.

###  6.2 Effective Forecast Horizon

Long forecast horizons in any statistical model can be challenging. This is
due to the recursive dependence that is common in time series. That is, models
typically will recycle outputs as future inputs in any out-of-sample
forecasts. As such, the error in any given point estimate will propagate to
future estimates, which limits the predictability of a given model. In the
context of chaotic time series, the Lyapunov time becomes a natural scale on
which to characterize an effective forecast horizon. That is, the time length
on which one can achieve sufficiently reliable forecasts. Typically, data sets
that exhibit more chaos in temporal dynamics will have a larger maximum
Lyapunov exponent, which effectively limits the real time range in predictive
power. Pathak et al. (2018) utilized the reservoir computing framework in
predicting spatio-temporally chaotic systems for the Kuramoto-Sivashinsky
equation. Their scheme achieves low prediction error for roughly 8 Lyapunov
times. As such, for comparative purposes we used $\Lambda t_{\text{max}}=8$ as
a benchmark. Our Opt-ESN forecasts achieved low NMSE beyond the 8 Lyapunov
time benchmark in each of the six simulation models (summarized in Table 2).
Beyond this point, it becomes more likely that the forecasts become
unreliable. In future work however, we may look into quantifying an asymptotic
upper limit on the model’s effective forecast horizon in terms of the Lyapunov
time.

##  7 Conclusion

We introduced the Opt-ESN model to forecast the mass accretion in protostellar
disk evolution. This model exploits the stochastic nature of echo state
networks and introduces its use in time-domain astronomy. We applied our model
to a series of synthetic mass accretion data sets simulated by solving the
hydrodynamical equations for a protostellar disk. The model achieved
predictions with a low normalized mean square error (NMSE) ($\sim 10^{-5}$ to
$10^{-3}$) for forecasts ranging between 0.099 to 3.793 kyr. Additionally, the
model successfully resolved the occurrence of an episodic burst with low NMSE
when we added 1% and 5% of noise to the data. However, our findings also
suggest that the model is not immune to degradation under scenarios of data
limitation.

Our implementation demonstrates the predictive capabilities of Opt-ESN when
applied to time series data. As we transition into a new era of time domain
astronomy, understanding and developing robust statistical time series models
is becoming increasingly important. Importantly, the scientific return of our
work goes much beyond its application to observations in the optical domain.
There may be synergies with observations made in the radio domain (e.g., fast
radio bursts) with facilities like the Canadian Hydrogen Intensity Mapping
Experiment (CHIME) and the Australian Square Kilometer Array Pathfinder
(ASKAP). Likewise in the gravitational-wave arena, Opt-ESN can play a crucial
role in detecting/denoising black hole and neutron star merger signals
observed from facilities like the Laser Interferometer Gravitational-Wave
Observatory (LIGO) and Virgo interferometer as well as the upcoming Kamioka
Gravitational Wave Detector (KAGRA) and LIGO India.

## acknowledgments

S.B. is supported by a Discovery Grant from NSERC. S.A. is supported by the
NASA Postdoctoral Program (NPP). S.A. acknowledges that a portion of this
research was carried out at the Jet Propulsion Laboratory, California
Institute of Technology, under a contract with the National Aeronautics and
Space Administration (80NM0018D0004). E.I.V. acknowledges support by the
Ministry of Science and Higher Education of the Russian Federation (State
assignment in the field of scientific activity 2023, GZ0110/23-10-IF).

Supplemental material is provided here for reference. This includes some
background discussion on stationarity, as well as additional figures and
tables that summarize particular sets of results.

###  .1 Stationarity

A process is considered to be stationary if the underlying distribution is
constant over time (Brockwell & Davis, 2002; Hamilton, 2020). This effectively
translates to having constant values of its first four moments. Alternatively,
we define a time series to be covariance-stationary if it is only constant in
its first and second moments. More formally, let ${X_{t}}$ be a time series
with $\mathbb{E}(X_{t}^{2})<\infty$. The first and second moments are the mean
and covariance functions, respectively:

$\displaystyle\mu_{X}(t)$ $\displaystyle=\mathbb{E}(X_{t}),$ (13)
$\displaystyle\Gamma_{X}(r,s)$
$\displaystyle=\text{Cov}(X_{r},X_{s}),\quad\forall r,s\,.$ (14)

Thus, the process ${X_{t}}$ is said to be covariance stationary if
$\mu_{X}(t)$ is independent of $t$ and $\Gamma_{X}(t+h,t)$ is independent of
$t$ for all $h$. Furthermore, the autocovariance function must be even and
nonnegative definite (Brockwell & Davis, 2002). That is, for a real-valued
vector A having components $a_{i}$, we have

$\displaystyle\sum_{i,j}a_{i}\Gamma_{X}(i-j)a_{j}\geq 0\,.$ (15)

Additionally, a stationary process will have the roots of its characteristic
equation lie inside the unit circle (Patterson, 2011). That is, if the
underlying process has a unit root $\geq 1$, then it is nonstationary. Assume
the variable ${X_{t}}$ can be written as a $p^{\text{th}}$ order
autoregressive process:

$\displaystyle
X_{t}=\alpha_{1}X_{t-1}+\alpha_{2}X_{t-2}+\dots+\alpha_{p}X_{t-p}+\epsilon_{t}\,,$
(16)

where the innovations $\epsilon_{t}$ are uncorrelated with mean zero and
constant variance. If the characteristic equation

$\displaystyle\lambda^{p}-\alpha_{1}\lambda^{p-1}-\alpha_{2}\lambda^{p-2}-\dots-\alpha_{p}=0$
(17)

has roots $\lambda\geq 1$ of multiplicity $m$, then the process is
nonstationary with integration order $m$, denoted $I(m)$. There are several
approaches to assessing whether a time series is stationary or not. The
augmented Dickey-Fuller (ADF) test is among a popular set of unit root tests
for time series data. It involves testing the null hypothesis that a unit root
is present in the underlying process, making the time series nonstationary
(Patterson, 2011). The ADF test assumes the underlying process can be modeled
by

$\displaystyle\Delta X_{t}=\alpha_{0}+\alpha_{1}t+\rho
X_{t-1}+\gamma_{1}\Delta X_{t-1}+\dots+\gamma_{p-1}\Delta
X_{p-t+1}+\epsilon_{t}\,.$ (18)

The process would have a unit root if $\rho=1$ or alternatively be considered
stationary if $\rho<1$. As such, the test is carried out under the null
hypothesis that $\rho=1$ against the alternative that $\rho<1$, where the test
statistic $\text{DF}_{\rho}$ is given as

$\displaystyle\text{DF}_{\rho}=\frac{\hat{\rho}}{\text{SE}(\hat{\rho})}\,.$
(19)

Here, $\text{SE}(\hat{\rho})$ is the standard error and the value
$\text{DF}_{\rho}$ is compared to the respective critical value in the Dickey-
Fuller distribution. In our implementation, we take stationarity to mean
covariance-stationary and utilize the ADF test at 5% significance with
$p=12\times(N_{t}/100)^{1/4}$ (default setting in Python), where $N_{t}$ is
the number of observations.

###  .2 Opt-ESN Model Outputs

Figure 12 demonstrates the Opt-ESN prediction of the fluctuating component of
the simulation data.

Figure 12: The Opt-ESN prediction on the fluctuating component of the
simulation data. The solid black line is representative of the simulation data
used in training/validation, the dotted black is the out-of-time testing
segment of the data and the red line is the network’s prediction.

Figure 13 demonstrates the Opt-ESN prediction of the deterministic component
of the simulation data.

Figure 13: The Opt-ESN prediction on the deterministic component of the
simulation data. The solid black line is representative of the simulation data
used in training/validation, the dotted black is the out-of-time testing
segment of the data and the red line is the network’s prediction.

###  .3 Performance Assessment

The tables below summarize the model performance metrics. Table 2 contains the
hyperparameters and forecast performances for each model component
(fluctuating, deterministic, and overall). Furthermore, Table 3 demonstrates
the performance metrics on forecasting the episodic burst under each of the
three training conditions.

Table 2: Opt-ESN hyperparameter settings and performance summary over each simulation model. Model | Component | Hyperparameters | $\Lambda$ | $\Lambda\cdot N_{t}$ | $t$ (kyr) | NMSE
---|---|---|---|---|---|---
$\varrho$ | $c_{r}$ | $\alpha$ | $\rho$ | Validation | Out-of-Time
| Fluctuating | 1.0000 | 0.1000 | 0.2000 | 0.5000 |  |  |  | $1.00\mathrm{e}{-4}$ | $4.05\mathrm{e}{-5}$
Model 26 | Deterministic | 0.1000 | 0.1000 | 0.2000 | 0.7000 | 0.208 | 10.41 | 1.978 | $1.33\mathrm{e}{-5}$ | $1.05\mathrm{e}{-5}$
| Overall |  |  |  |  |  |  |  | $5.11\mathrm{e}{-5}$ | $1.86\mathrm{e}{-5}$
| Fluctuating | 1.0000 | 0.1000 | 0.9000 | 0.9900 |  |  |  | $3.75\mathrm{e}{-4}$ | $2.91\mathrm{e}{-4}$
Model 27 | Deterministic | 0.1886 | 0.1000 | 0.4655 | 0.1957 | 0.201 | 10.07 | 0.975 | $3.15\mathrm{e}{-4}$ | $1.01\mathrm{e}{-3}$
| Overall |  |  |  |  |  |  |  | $3.84\mathrm{e}{-4}$ | $8.15\mathrm{e}{-4}$
| Fluctuating | 0.0769 | 0.0100 | 0.1117 | 0.6774 |  |  |  | $1.68\mathrm{e}{-3}$ | $2.32\mathrm{e}{-3}$
Model 29 | Deterministic | 0.5000 | 0.1000 | 0.2000 | 0.5000 | 0.184 | 9.22 | 0.341 | $1.99\mathrm{e}{-4}$ | $2.94\mathrm{e}{-4}$
| Overall |  |  |  |  |  |  |  | $1.34\mathrm{e}{-3}$ | $2.28\mathrm{e}{-3}$
| Fluctuating | 0.0298 | 0.0100 | 0.2680 | 0.3306 |  |  |  | $1.78\mathrm{e}{-4}$ | $7.97\mathrm{e}{-5}$
Model 30 | Deterministic | 0.0783 | 0.0100 | 0.1970 | 0.2960 | 0.162 | 8.12 | 3.793 | $1.53\mathrm{e}{-5}$ | $2.15\mathrm{e}{-5}$
| Overall |  |  |  |  |  |  |  | $1.83\mathrm{e}{-4}$ | $5.93\mathrm{e}{-5}$
| Fluctuating | 0.5000 | 0.1000 | 0.5000 | 0.8500 |  |  |  | $1.06\mathrm{e}{-3}$ | $1.23\mathrm{e}{-3}$
Model 31 | Deterministic | 0.5000 | 0.1000 | 0.5000 | 0.2000 | 0.199 | 9.94 | 0.312 | $1.58\mathrm{e}{-4}$ | $8.95\mathrm{e}{-4}$
| Overall |  |  |  |  |  |  |  | $1.15\mathrm{e}{-3}$ | $1.91\mathrm{e}{-3}$
| Fluctuating | 0.0673 | 0.0100 | 0.3224 | 0.2933 |  |  |  | $7.64\mathrm{e}{-4}$ | $9.50\mathrm{e}{-4}$
Model 32 | Deterministic | 0.0481 | 0.0100 | 0.1590 | 0.3532 | 0.192 | 9.61 | 0.099 | $1.02\mathrm{e}{-4}$ | $9.98\mathrm{e}{-4}$
| Overall |  |  |  |  |  |  |  | $2.85\mathrm{e}{-4}$ | $1.18\mathrm{e}{-3}$
Table 3: Opt-ESN model settings used in episodic burst predictions. Condition | Component | Hyperparameters | NMSE | Training length | Level of noise
---|---|---|---|---|---
$\varrho$ | $c_{r}$ | $\alpha$ | $\rho$
| Fluctuating | $0.6661$ | $0.1000$ | $0.1836$ | $0.6394$ |  $1.64$e$-2$ |  |
(1) | Deterministic | $0.8185$ | $0.1000$ | $0.2583$ | $0.2953$ |  $2.46$e$-4$ | 10000 | 0%
| Overall |  |  |  |  |  $8.29$e$-3$ |  |
| Fluctuating | 0.0390 | 0.1000 | 0.1535 | 0.4175 |  $2.24$e$-2$ |  |
| Deterministic | 0.2091 | 0.1000 | 0.3273 | 0.5917 |  $1.81$e$-3$ | 5000 | 0%
| Overall |  |  |  |  |  $1.07$e$-2$ |  |
| Fluctuating | 0.0474 | 0.1000 | 0.4446 | 0.1319 |  $5.16$e$-2$ |  |
(2) | Deterministic | 0.6935 | 0.1000 | 0.7389 | 0.3642 |  $3.10$e$-3$ | 2500 | 0%
| Overall |  |  |  |  |  $2.25$e$-2$ |  |
| Fluctuating | 0.7489 | 0.1000 | 0.0102 | 0.9848 |  $1.14$e$+5$ |  |
| Deterministic | 0.0184 | 0.1000 | 0.4612 | 0.1193 |  $1.74$e$+0$ | 1000 | 0%
| Overall |  |  |  |  |  $6.38$e$+4$ |  |
| Fluctuating | $0.3026$ | $0.1000$ | $0.1755$ | $0.8945$ |  $2.58$e$-2$ |  |
| Deterministic | $0.6681$ | $0.1000$ | $0.4042$ | $0.2133$ |  $1.27$e$-3$ | 10000 | 1%
(3) | Overall |  |  |  |  |  $1.84$e$-2$ |  |
| Fluctuating | $1.0000$ | $0.1000$ | $0.3000$ | $0.9000$ |  $6.81$e$-2$ |  |
| Deterministic | $0.7527$ | $0.1000$ | $0.0296$ | $0.8296$ |  $1.36$e$-1$ | 10000 | 5%
| Overall |  |  |  |  |  $1.47$e$-1$ |  |

## References

  * Audard et al. (2014) Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 387, doi: http://doi.org/10.2458/azu_uapress_9780816531240-ch01710.2458/azu_uapress_9780816531240-ch017
  * Auddy et al. (2022) Auddy, S., Dey, R., Lin, M.-K., Carrera, D., & Simon, J. B. 2022, ApJ, 936, 93, doi: http://doi.org/10.3847/1538-4357/ac7a3c10.3847/1538-4357/ac7a3c
  * Auddy et al. (2021) Auddy, S., Dey, R., Lin, M.-K., & Hall, C. 2021, The Astrophysical Journal, 920, 3, doi: http://doi.org/10.3847/1538-4357/ac151810.3847/1538-4357/ac1518
  * Auddy & Lin (2020) Auddy, S., & Lin, M.-K. 2020, ApJ, 900, 62, doi: http://doi.org/10.3847/1538-4357/aba95d10.3847/1538-4357/aba95d
  * Bai (2013) Bai, X.-N. 2013, ApJ, 772, 96, doi: http://doi.org/10.1088/0004-637X/772/2/9610.1088/0004-637X/772/2/96
  * Balbus (2003) Balbus, S. A. 2003, ARA&A, 41, 555, doi: http://doi.org/10.1146/annurev.astro.41.081401.15520710.1146/annurev.astro.41.081401.155207
  * Bloom & Richards (2012) Bloom, J. S., & Richards, J. W. 2012, in Advances in Machine Learning and Data Mining for Astronomy, ed. M. J. Way, J. D. Scargle, K. M. Ali, & A. N. Srivastava, 89–112, doi: http://doi.org/10.48550/arXiv.1104.314210.48550/arXiv.1104.3142
  * Brockwell & Davis (2002) Brockwell, P. J., & Davis, R. A. 2002, Introduction to time series and forecasting (Springer)
  * Caratti o Garatti et al. (2017) Caratti o Garatti, A., Stecklum, B., Garcia Lopez, R., et al. 2017, Nature Physics, 13, 276, doi: http://doi.org/10.1038/nphys394210.1038/nphys3942
  * Contreras Peña et al. (2019) Contreras Peña, C., Naylor, T., & Morrell, S. 2019, MNRAS, 486, 4590, doi: http://doi.org/10.1093/mnras/stz101910.1093/mnras/stz1019
  * Cybenko (1989) Cybenko, G. 1989, Mathematics of control, signals and systems, 2, 303
  * Das & Basu (2022) Das, I., & Basu, S. 2022, MNRAS, 514, 5659, doi: http://doi.org/10.1093/mnras/stac165410.1093/mnras/stac1654
  * Eckmann et al. (1986) Eckmann, J.-P., Kamphorst, S. O., Ruelle, D., & Ciliberto, S. 1986, Physical Review A, 34, 4971
  * Elbakyan et al. (2016) Elbakyan, V. G., Vorobyov, E. I., & Glebova, G. M. 2016, Astronomy Reports, 60, 879, doi: http://doi.org/10.1134/S106377291610001210.1134/S1063772916100012
  * Hamilton (2020) Hamilton, J. D. 2020, Time series analysis (Princeton university press)
  * Jaeger (2007) Jaeger, H. 2007, scholarpedia, 2, 2330
  * Kim & King (2020) Kim, T., & King, B. R. 2020, Neural Computing and Applications, 32, 17769
  * Liu et al. (2018) Liu, X., Chen, M., Yin, C., & Saad, W. 2018, in 2018 17th IEEE International Conference on Machine Learning and Applications (ICMLA), IEEE, 443–448
  * Lukoševičius (2012) Lukoševičius, M. 2012, Neural Networks: Tricks of the Trade: Second Edition, 659
  * Machida & Basu (2019) Machida, M. N., & Basu, S. 2019, ApJ, 876, 149, doi: http://doi.org/10.3847/1538-4357/ab18a710.3847/1538-4357/ab18a7
  * Meyer et al. (2017) Meyer, D. M. A., Vorobyov, E. I., Kuiper, R., & Kley, W. 2017, MNRAS, 464, L90, doi: http://doi.org/10.1093/mnrasl/slw18710.1093/mnrasl/slw187
  * Pathak et al. (2018) Pathak, J., Hunt, B., Girvan, M., Lu, Z., & Ott, E. 2018, Physical review letters, 120, 024102
  * Patterson (2011) Patterson, K. 2011, Unit root tests in time series volume 1: key concepts and problems (Springer)
  * Rocha-Solache et al. (2022) Rocha-Solache, A., Rodríguez-Montoya, I., Sánchez-Argüelles, D., & Aretxaga, I. 2022, ApJS, 260, 15, doi: http://doi.org/10.3847/1538-4365/ac525910.3847/1538-4365/ac5259
  * Salehinejad et al. (2017) Salehinejad, H., Sankar, S., Barfett, J., Colak, E., & Valaee, S. 2017, arXiv e-prints, arXiv:1801.01078, doi: http://doi.org/10.48550/arXiv.1801.0107810.48550/arXiv.1801.01078
  * Schäfer & Zimmermann (2006) Schäfer, A. M., & Zimmermann, H. G. 2006, in International Conference on Artificial Neural Networks, Springer, 632–640
  * Suzuki et al. (2016) Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74, doi: http://doi.org/10.1051/0004-6361/20162895510.1051/0004-6361/201628955
  * Vorobyov & Basu (2005) Vorobyov, E. I., & Basu, S. 2005, ApJ, 633, L137, doi: http://doi.org/10.1086/49830310.1086/498303
  * Vorobyov & Basu (2006) —. 2006, ApJ, 650, 956, doi: http://doi.org/10.1086/50732010.1086/507320
  * Vorobyov & Basu (2007) —. 2007, MNRAS, 381, 1009, doi: http://doi.org/10.1111/j.1365-2966.2007.12321.x10.1111/j.1365-2966.2007.12321.x
  * Vorobyov & Basu (2009) —. 2009, MNRAS, 393, 822, doi: http://doi.org/10.1111/j.1365-2966.2008.14376.x10.1111/j.1365-2966.2008.14376.x
  * Vorobyov & Basu (2010) —. 2010, ApJ, 719, 1896, doi: http://doi.org/10.1088/0004-637X/719/2/189610.1088/0004-637X/719/2/1896
  * Vorobyov & Basu (2015) —. 2015, ApJ, 805, 115, doi: http://doi.org/10.1088/0004-637X/805/2/11510.1088/0004-637X/805/2/115
  * Vorobyov et al. (2017) Vorobyov, E. I., Elbakyan, V., Hosokawa, T., et al. 2017, A&A, 605, A77, doi: http://doi.org/10.1051/0004-6361/20163035610.1051/0004-6361/201630356
  * Vorobyov et al. (2020) Vorobyov, E. I., Khaibrakhmanov, S., Basu, S., & Audard, M. 2020, A&A, 644, A74, doi: http://doi.org/10.1051/0004-6361/20203908110.1051/0004-6361/202039081
  * Vulpiani et al. (2009) Vulpiani, A., Cecconi, F., & Cencini, M. 2009, Chaos: from simple models to complex systems, Vol. 17 (World Scientific)