* [Sch91] J. Schaeffer.  Global existence of smooth solutions to the Vlasov–Poisson system in three dimensions.  Comm. Partial Differential Equations, 16(8-9):1313–1335, 1991.
  * [Sch04] J. Schaeffer.  Steady states in galactic dynamics.  Arch. Ration. Mech. Anal., 172(1):1–19, 2004.
  * [Smu16] J. Smulevici.  Small data solutions of the Vlasov–Poisson system and the vector field method.  Ann. PDE, 2(2):Art. 11, 55, 2016.
  * [SS06] Ó. Sánchez and J. Soler.  Orbital stability for polytropic galaxies.  Ann. Inst. H. Poincaré C Anal. Non Linéaire, 23(6):781–802, 2006.
  * [Vla68] A. A. Vlasov.  The vibrational properties of an electron gas.  Soviet Physics Uspekhi, 10(6):721–733, 1968.
  * [Wan22] X. Wang.  Propagation of regularity and long time behavior of the 3$D$ massive relativistic transport equation II: Vlasov–Maxwell system.  Comm. Math. Phys., 389(2):715–812, 2022.
  * [WZ11] J. Wunsch and M. Zworski.  Resolvent estimates for normally hyperbolic trapped sets.  Ann. Henri Poincaré, 12(7):1349–1385, 2011.


# Employing Typicality in Optimal Control Theory

Aviv Aroch  The Institute of Chemistry, The Hebrew University of Jerusalem,
Jerusalem 9190401, Israel  Shimshon Kallush  Ronnie Kosloff  The Institute of
Chemistry, The Hebrew University of Jerusalem, Jerusalem 9190401, Israel

###### Abstract

Controlling the dynamics of quantum systems is a crucial task in quantum
science and technology. Obtaining the driving field that transforms the
quantum systems to its objective is a typical control task. This task is hard,
scaling unfavorably with the size of Hilbert space. To tackle this issue we
employ typicality to assist in finding the control field for such systems. To
demonstrate the method we choose the control task of cooling the fine
structure states of the AlF molecule, from relatively high temperatures which
results in large Hilbert space. Using quantum typicality, we demonstrate that
we can simulate an ensemble of states, enabling a control task addressing
simultaneously many states. We employ this method to find a control field for
cooling molecules with large number of internal sates, corresponding to high
initial temperatures.

##  I Introduction

Controlling quantum phenomena has been a goal in quantum physics and
chemistry. Quantum control theory addresses this topic and has evolved rapidly
over the last three decades [1]. One of the main goals of quantum control
theory is to establish a series of systematic methods to manipulate and
control quantum systems. Quantum control theory has been implemented in
physical chemistry, atomic and molecular physics, and quantum optics. In much
of quantum control theory, the controllability of quantum systems is a
fundamental issue [2]. Controllability concerns the existence of a control
solution for a specific task. This problem has practical importance since it
closely connects with the universality of quantum computation and the
possibility of achieving atomic or molecular scale transformations. For
finite-dimensional quantum systems, the controllability criteria may be
expressed in terms of the structure and rank of corresponding Lie groups and
Lie algebras [3, 4]. This method allows for a mathematical treatment of the
problem in the case of closed quantum systems.

The existence of a controllable task does not hint at how to obtain the
control field. Iterative schemes have been developed based on constrained
optimization has been developed [5, 6, 7, 8, 9, 10, 11]. Controlling systems
with Large Hilbert space is a challenging problem mainly due to its
computational complexity.

The computational effort measures the complexity of the problem we are trying
to solve. In our case, the problem is twofold; our primary step requires
solving the time-dependent Schrödinger equation for each state. The complexity
scales between $M\times N\log N$ and $M\times N^{2}$ where $N$ is the size of
Hilbert space of our system and $M$ the number of time-steps which in turn
scale as $M\propto O(\Delta E\times t)$ where $t$ is the time interval and
$\Delta E$ the energy range [12]. The number of iterations for solving a
unitary control problem scales like $K!$ where $K$ is the size of the
transformation [7]. The complexity class of the cooling transformation, which
addresses $K$ states simultaneously, is not known and could be between
polynomial in $K$ and factorial $K!$. Altogether the control problem is
computationally highly complex [13].

This issue requires establishing methods that reduce the number of states $K$
needed for optimization and the complexity. To achieve this task, we will
employ the properties of quantum typicality.

Quantum typicality states that a single quantum state can typically well
describe local expectation values of a quantum system. This statement applies
to Schrödinger type dynamics in high dimensional Hilbert spaces. As a
consequence, individual dynamics of expectation values converge to the
ensemble’s average [14]. We will employ quantum typicality by using random
phase wavefunctions for control problems.

The random phase wave function method uses an ensemble of pure states, which
creates an efficient representation of the mixed state of the entire system.
The convergence of the RPWF method becomes faster as the size of the Hilbert
space (and the number of random phases) increases, namely with the rise of
initial temperature.

The present study aims to develop an optimal control algorithm to cool the
molecular internal degrees of freedom, a multi-state control task. The control
objective is to increase the system’s purity ${\cal
P}=Tr\\{\hat{\rho}^{2}\\}$. Alternatively, cooling can be defined as lowering
the effective temperature $T_{eff}$ defined by the von Neumann entropy of the
ensemble [15]. $T_{eff}$ is the temperature of a passive Gibbs state with the
same von Neumann entropy [16].

The field of cold and ultra-cold molecules has rapidly grown in the last
decade. Cold molecules have an essential role in many active areas in science
amongst tests of fundamental physics, cold chemistry, quantum technologies
(QT), and quantum information (QI)[17, 18, 19, 20, 21].

There are many ways of cooling particles; one of the most effective is
radiative cooling. The basic technique is to locate a closed loop of
stimulated excitation and spontaneous emission. The cooling is achieved by
entropy removal by spontaneous emission. This technique has been employed to
reduce the temperature of translational as well as internal degrees of freedom
of atomic and molecular species [22, 23, 24].

A crucial consideration in cooling internal degrees of freedom in molecules is
enforcing a close transition cycle. Excitation may result in the molecule
ending in a different state outside the closed-loop after decay; thus, cooling
molecules by laser excitation is complicated, and only a few examples exist.
However, the molecular complexity makes cold molecules applicable for a
broader range than atoms.

We choose optimal control theory to overcome the molecular complexity and
enforce closed-loop solutions. Controlled laser fields are employed to remove
frequencies that damage the required transition. Control can steer a quantum
system from its initial state to its final one (target state). Optimal Control
Theory (OCT) was created to do this task with maximum fidelity (defined by the
user).

To control a specific system, we first need to describe the state and its
evolution in time. In our system, we can define the rotational state employing
the eigenstates of a rigid rotor as a complete basis. The state is propagated
as an isolated quantum system evolving by a unitary transformation.

Optimal Control Theory (OCT) is employed for an isolated system to obtain the
field leading to the target state. However, this alone will not lead us into a
colder state (pure state); entropy is invariant under unitary transformation.
Therefore cooling requires dissipation that can change the entropy of our
quantum system.

Assuming that the target of the process is to reach a pure, single state, any
proposed mechanism for the process has to maintain the population of the
single-target state while allowing the population of all other states to
repopulate selectively. However, as was shown by [25, 26], the speciﬁc choice
of the precooling transformation is a subtle issue. As implied by the ergodic
theorem [27], for any initial state under multiple cycles of a given ﬁeld-
driven unitary transformation and subsequent decay, the ﬁnal state will be the
invariant under the whole transformation. After many excitation-relaxation
cycles, the memory of the initial state will be erased. The entire
transformation can be described by the following, let ${\cal U}$ be the
unitary super operator and ${\cal D}$ be dissipative super operator; then
there is a stationary state $\hat{\rho}_{SS}$ that will obey

$\lim_{n\rightarrow\infty}\left({\cal U}{\cal
D}\right)^{n}\hat{\rho}=\hat{\rho}_{ss}$ (1)

where $\hat{\rho}$ can be any initial state, and ${\cal
U}\bullet=\hat{U}\bullet\hat{U}^{\dagger}$. That is, under a given ${\cal U}$
and ${\cal D}$, the system will ﬁnally evolve from any state into the single
stationary $\hat{\rho}_{ss}$. In this way, we look for a unitary control
transformation that under a known dissipator will lead into a purer state
[22].

The ultimate goal of molecular cooling is to reach a large sample of pure
state molecules that can be transformed to a BEC. This task can only be
achieved in stages. The initial sample should be as pure as possible, pure
from other isotopes, other quantum states (hotter ones in our case), etc., is
crucial. A crucial intermediate step is to obtain a molecular Magnetic Optical
Trap (MOT). For achieving this step, an ensemble of a single rovibrational
state should be created.

Our control objective is designed to enhance the initial preparation; it can
be used as a purification method to add molecules that are not resonant with
MOT transition to be added and thus increase the number density of the initial
ensemble.

Specifically, we choose to cool the AlF molecule rotationally. Spectroscopic
measurements and detailed analysis [28], have shown that such a task is
feasible using optimal control theory. This molecule belongs to a family of
which the vibrational manifolds are closed due to a Frank-Condon coefficient
close to one. Our obstacle is the large initial number of rotational states.

The present theoretical study aims to obtain a high fidelity control task with
low computational effort. We have implemented this method on the fine-
structure levels of the AlF molecule with a large total orbital angular
momentum at 30K. We will employ optimal control theory to ﬁnd broadband-shaped
pulses, steering the system into a colder state.

The paper is arranged as follows: Section II describes the model in which we
implemented our tools, in section III describes the theoretical tools used for
achieving control. Section IV presents the results, which are discussed and
summarized in the concluding Sec. V.

##  II The Model

The model we employ is the rovibrational structure of the AlF molecule. The
state of the systems is defined by the combined density operator of two
electronic surfaces:

$\\\
{\hat{\rho}}={\hat{\rho}_{g}}\otimes{\hat{P}_{g}}+{\hat{\rho}_{e}}\otimes{\hat{P}_{e}}+{\hat{\rho}_{c}}\otimes{\hat{S}_{+}}+\hat{\rho}_{c}^{*}\otimes{\hat{S}_{-}}\
$ (2)

Where $\hat{P}_{g/e}$ are the projection operators of the ground and excited
electronic state, $\hat{S}_{+/-}$ are the electronic raising and lowering
operators, $\rho_{g/e}$ are the density operator for the rovibrational
ensemble within the ground and excited electronic states, and $\hat{\rho}_{c}$
is the density operator of the nuclear coherence between the surfaces.

For AlF, we assume an initial temperature of $\sim 30^{0}K$, for which the the
population of the states occupies up to $J=11$ (144 sub-levels).

The Liouville von Neumann equation governs the evolution of the system:

$\frac{{d\hat{\rho}}}{{dt}}=-\frac{i}{\hbar}[\hat{H},\hat{\rho}]+{\cal
L}_{D}(\hat{\rho})$ (3)

The first term is the coherent dynamical part governed by the Hamiltonian and
the second is the dissipative part of the dynamics describing spontaneous
emission. This equation represents the dynamics of an open quantum system.

For optical transitions with multiple pulses, there is a distinct timescale
separation between the light-induced step, which occurs in less than a
picosecond and is unitary, the incoherent decay which occurs in tens of
nanoseconds, and the pulse repetition rate, which is in the $MHz$ to $KHz$.
Each cooling cycle could be separated within this picture into two parts: (1)
The short-time interaction of the external field and the molecular system.
Since this step is unitary, the density operator can be decomposed to energy
eigenstates, and each component can be computed in a wave function framework.
Then (2) a slow and field-free, spontaneous decay takes place. In this step,
the coherences developed between energy eigenstates during the laser-
controlled stage are erased.

The Hamiltonian, which governs the unitary part of the dynamics, can be
written as :

$\begin{array}[]{l}{\hat{H}_{t}={{\hat{H}}_{0}}}+{{\hat{V}}_{t}}\\\
{{\hat{H}}_{0}}={{\hat{H}}_{g}}\otimes{{\hat{P}}_{g}}+{{\hat{H}}_{e}}\otimes{{\hat{P}}_{e}}.\end{array}\
$ (4)

Where ${\hat{H}}_{g/e}$ is the ground and excited rotatioanl Hamiltonian. The
interaction of the system with light, assumed to be linearly polarized to the
lab $z$ axis is described by ${\hat{V}}_{t}$:

${{\hat{V}}_{t}}=-\hat{\mu}_{z}\otimes({{\hat{S}}_{+}}\varepsilon_{z}(t)+{{\hat{S}}_{-}}{\varepsilon_{z}(t)^{*}})$
(5)

Where $\hat{\mu}_{z}$ is the transition dipole moment along the $z$ spatial
direction and $\varepsilon_{z}(t)$ represents the time-dependent field along
the same direction.

For mildly cold temperatures ($T_{initial}\leq 30K$) and high vibrational
frequency, we can assume that the molecules are initially in their ground
$v=0$ state. ${\hat{H}_{g/e}}$ are the field-free rotational Hamiltonians for
the ground vibrational state. Moreover, for the $AlF$ and chemically similar
molecules, vibrational excitations in the electronic transition are negligible
due to the highly restricting FC factors [29]. As a result the modeling is
thus restricted to $v=0$.

Under the Hund’s case $a$, applicable to our case, the rotational states of
the model are expanded by the symmetric top basis [30]

$|J,\Omega
M\rangle=\left[\frac{2J+1}{4\pi}\right]^{\frac{1}{2}}D_{M,\Omega}^{J}(\phi,\theta,0)$
(6)

Where $J$ is the total molecular angular momentum, $M$ and $\Omega$ are
projections on the spatial (Z) and molecular (z) axes, respectively. Here,
$D_{M\Omega}^{J}$ is the rotational tensor. For the ground electronic
${}^{1}\Sigma$ state, the projection of the spatial electronic angular
momentum on the molecular axis is 0 (L=S=0), therefor $\Lambda=0$, and J=N=R.
The rotational Hamiltonian becomes

$\hat{H}_{\text{rot}}(r)=B\hat{\bm{R}}^{2}=B\hat{\bm{J}}^{2}$ (7)

where $B$ is a rotational constant, and $\hat{R}$ is the nuclear rotational
angular momentum operator, which is equal to ${\bm{R}}={J}-{L}-{S}$ , L is the
electronic orbital angular momentum and S is the electronic spin angular
momentum.

${}^{1}\Sigma$ energies are then

$E(^{1}{\Sigma};J)={B}J(J+1)$ (8)

and its corresponding eigenstates can be defined by $|J\Omega M\rangle$ when
$\Omega=0$.
In the excited ${}^{1}\Pi$ state $\Lambda=\pm 1$ S=0, and thus $\Omega=\pm 1$.

The rotational Hamiltonian for the excited state, ${}^{1}\Pi$, is

$\hat{H}_{\text{rot}}(r)=B\hat{\bm{R}}^{2}=B(\hat{\bm{J}}-\hat{\bm{L}})^{2}$
(9)

its energies are

$E(^{1}{\Sigma};J)={B}[J(J+1)-1]$ (10)

the corresponding eigenstates are $|J,\Omega M\rangle$ where $\Omega=\pm{1}$.

The transitions between the two electronic states
${}^{1}\Sigma\rightarrow{}^{1}\Pi$ are dictated by dipole selection rules
$(\Delta J=0,\pm 1)$, denoted as R,Q,and P branches,respectively.The coupling
elements can be found by calculating the overlap of any two eigenvectors with
the dipole operator,

$\displaystyle\int
D_{m^{\prime}\Omega^{\prime}}^{j^{\prime}}(\theta,\phi,0)D_{0q}^{1}(\theta,\phi,0)D_{m\Omega}^{j}(\theta,\phi,0)d\Omega$
(11) $\displaystyle=8\pi\left(\begin{array}[]{ccc}j&1&j^{\prime}\\\
-m&0&m^{\prime}\end{array}\right)\left(\begin{array}[]{ccc}j&1&j^{\prime}\\\
-\Omega&q&\Omega^{\prime}\end{array}\right)$

where $\mu_{0q}$ is the transition dipole moment proportional to
$D_{0q}^{1}(\theta,\phi,0)$. The value of q is determined for a given
transition case and is equal to $q=\Omega-\Omega^{\prime}$.

At thermal equilibrium, the initial state $\rho_{eq}$ is characterized by
thermally distributed in the quantum rotational states

$\hat{\rho}_{eq}=\frac{1}{Z}\sum_{j}e^{-\beta E_{g}^{j}}|J0M\rangle\langle
J0M|$ (12)

where $\beta=1/k_{b}T$ and the sets $|J\Omega M\rangle$ and
$\left\\{E_{j}\right\\}$ are the eigenstates and eigenenergies of the system,
and Z is the partition function.

The dissipating part of the dynamics is generated by the Liovillian super
operator ${\cal L}_{D}$, Eq. (3). Integrating in time leads to the transition
map $\Lambda_{t}=e^{{\cal L}_{D}t}$. Assuming that the timescale between
pulses is longer than the spontaneous emission, the transition map ${\cal D}$
of the spontaneous emission can be defined. The matrix elements of ${\cal D}$
are from a given excited state energy eigenstate to the ground state manifold
of states. They can be calculated employing Fermi’s golden rule:

$\Gamma_{i\rightarrow f}=\frac{2\pi}{\hbar}\left|\left\langle
f\left|\mu\right|i\right\rangle\right|^{2}\rho\left(E_{f}\right)$ (13)

where $\left\langle f\left|H^{\prime}\right|i\right\rangle$ is the matrix
element of the electronic transition dipole between the final and initial
states, $\rho\left(E_{f}\right)$ is the density of states in vacuum (number of
continuum states divided by $dE$ in the infinitesimally small energy interval
$E$ to $E+dE$ ) at the energy $E_{f}$ of the final states. With Eq. (13) we
have defined the decay rate ($\sim 10^{-6}~{}second$), it is important to note
that because the rotattional selection rules ,in our system, dictates narrow
band transitions (all transitions are almost equal), $\omega^{3}$ can be
neglected. Assuming that we wait sufficient time and normalize it; we get the
decay probability of our system.

It is important to note that the coherent step yields only unitary
transformations and therefore does not change the purity when the latter decay
step is known for the loss of purity and thus changes the system’s
temperature.

We can use the defined Full map transformation Eq. (1) and get a final state
that is invariant under the whole transformation. After many excitation-
relaxation cycles, the memory of the initial state will be erased and the
system will finally evolve into the single stationary $\rho_{ss}$.

One can obtain the state $\rho_{ss}$ by diagonalization of Eq. (1). The
eigenstate then gives the stationary state with a unit eigenvalue, while the
next eigenvalue indicates the system’s convergence rate to the steady-state.

Using the ergodic theory compels the system to be closed. Since ${\cal D}$ is
fixed, the task is to find the unitary transformation ${\cal U}$ that is will
lead to the designed stationary state.

To associate an effective temperature to the state, we employ the von-Neuman
entropy to scale the purity and define the effective temperature. The idea
comes from information theory, where the entropy is related to the probability
distribution of an ensemble [15]. The entropy is defined as:

${S_{VN}}=-tr\\{\hat{\rho}\ln\hat{\rho}\\}\leq-\sum\limits_{j}{{P_{j}}}ln{P_{j}}$
(14)

where $\hat{\rho}$ is the system’s density matrix, and $P_{j}$ is the
probability to be in the energy eigenstate $j$. This is the only contribution
to the entropy, assuming that quantum coherences do not survive the
spontaneous emission incoherent step. Equality will be obtained when the
system is diagonalized in the energy domain. It is important to note that
entropy is invariant under unitary transformation. Therefore any steady-state
reached after cooling can be transformed by unitary transformation to a
passive state with the same entropy [31]. To define a temperature of any non-
thermal state, we associate it with the temperature of a thermal state with
the same vN entropy [16].

##  III Methods

###  III.1 Optimal Control Theory (OCT)

Quantum Optimal Control Theory (OCT) is a branch of coherent control, a
quantum mechanical based method for controlling dynamical processes. The basic
principle is to control quantum interference phenomena typically by shaping
the phase of laser pulses [32, 33, 5, 6]. OCT is formulated as a maximization
problem, and seeks a time dependent field that maximizes the expectation value
of an operator in final time.

Consider a quantum system in an initial state:
${\hat{\rho}_{0}}=\sum\limits_{k=1}^{N}p_{k}|\psi_{k}^{0}\rangle\langle\psi_{k}^{0}|$,
where the set $\\{\psi_{k}^{0}\\}$ is energy eigenstates of the system, and
$N$ the size of Hilbert space. The control will seek a field that maximizes
the expectation value of the operator $\hat{O}$ at final time T:

${J_{\max}}(\varepsilon)\equiv\sum\limits_{k=1}^{N}{p_{k}\langle{\psi_{k}}(T)|\hat{O}\left|{{\psi_{k}}(T)}\right\rangle}$
(15)

where ${\Psi_{i}}(T)$ describes the state that results from the interaction of
the system with the field $\varepsilon$at the final time $T$. The governing of
the dynamics of the system by the Schrödinger equation
$i\frac{\partial}{\partial t}|\psi\rangle={\hat{H}}|\psi\rangle$. The quantum
dynamics is enforced by adding an additional cost term to the functional,
according to the Lagrange-multiplier method:

${J_{con}}=\sum\limits_{k=1}^{N}{-2{\rm{Re}}\int\limits_{0}^{T}{\langle{\chi_{k}}(t)|\frac{d}{{dt}}+i\hat{H}(t)}\left|{{\psi_{k}}(t)}\right\rangle
dt}$ (16)

where $\\{\left\langle{\chi_{k}}(t)\right|\\}$ are the set of time dependent
Lagrange function multipliers. To regularize the solution with a limitation
over the intensity another penalty term to the functional is added [34]:

${J_{penal}}(\varepsilon)=-{\alpha}{\int_{0}^{T}{\varepsilon^{2}(t)dt}}$ (17)

where $\alpha$ is a scalar Lagrange multiplier.

The overall object of maximization is the following functional:

$J=J_{max}+J_{penal}+J_{con}$ (18)

The maximization of the fitness $J$ is the control task. Functional
derivatives with respect to the various field are then taken resulting in the
following system of equations:

Each of the set of the $|\chi_{k}(t)\rangle$ Lagrange multipliers will obey a
time reversed Schrödinger equation:

$\frac{d{\left\langle{\chi_{k}}(t)\right|}}{dt}=i\left\langle{{\chi_{k}}(t)}\right|\hat{H}(t)$
(19)

with the boundary conditions: $|\chi_{k}(T)\rangle=\hat{O}|\Psi_{k}(T)\rangle$
[35].

The Krotov iterative method is applied to obtain a monotonic growth of of the
fitness $J$ at each iteration with the updated field so that:

$\varepsilon^{l+1}(t)=\varepsilon^{l(t)}-\frac{1}{\alpha}{{\sum\limits_{k=1}^{N}{{\mathop{\rm
Im}\nolimits}\left\langle{\chi_{k}^{l}(t)}\right|\hat{\mu}}\left|{\psi_{k}^{l+1}(t)}\right\rangle}}$
(20)

where $\bullet^{{}^{l}}(t)$ denote the quantity after the the field after the
$l-th$ iteration.

Note that the scheme of Eqs. (19) (20) is somewhat similar to the simultaneous
optimization scheme that is required for unitary transformations and quantum
gates. However, for cooling each cycle at the final transformation erases the
relative quantum phase between the various optimized set initial states. This
leaves the resulting fitness measurement at the level of classical transition
probability between the initial and final state and removes the need to
evaluate quantum phases.

For large systems, the convergence of OCT is fast in the first iterations but
later saturates or becomes stagnant. Therefore, accurate solutions for large
systems of the control problem are difficult to reach. A good initial guess
for the control field can speed convergence considerably. Ideally, a control
field of an easy-to-solve small quantum system could serve as an initial pilot
guess for large quantum systems [19].

The first step is to find the pilot field for which we seek a solution to the
control problem of dimention $K$ for a small number of random phase states
$L$. We will then explore the universality of the control field for the
increasing number of states. We will use optimal control theory, based on
Krotov’s algorithm, which guarantees monotonic convergence [36], to do so.

The difficulty in converging a control field to generate state-to-state
transitions can be related to the algebraic structure of the control
Hamiltonian. We have seen that when the initial or target states are
superpositions of generalized coherent state, no relation exists between the
control fields of such targets for different Hilbert space sizes.

###  III.2 Quantum Typicality

Typicality describes a property of a system where a typical state can present
an assembly of similar states. This set of states should have a narrow
distribution of some feature (e.g., drawn according to the same distribution,
sharing the same energy, etc.) and therefore yield a very limited distribution
of expectation values. The typical state will fit the expectation value of the
complete set of states.

Quantum typicality was first noted by Schrödinger and von Neumann when they
were trying to incorporate statistical mechanics with quantum mechanics. They
inferred that the wavefunction of a complex system can have statistical
properties.

In their approach, when discussing thermalization in isolated quantum systems,
one should focus on physical observables instead of wave functions or density
matrices describing the entire system. This approach is similar to Eigenvalue
Thermaliziation Hypothesis (ETH), in which the focus is put on macroscopic
observables and “typical” configuration.

Eigenvalue Thermaliziation Hypothesis (ETH) implies that the expectation
values of local observables and their fluctuations in isolated quantum systems
relax to (nearly) time-independent values that can be described using
traditional statistical mechanics ensembles. This has been verified in several
quantum lattice systems and, according to ETH, should occur in generic many-
body quantum systems. ETH states that the eigenstates of generic quantum
Hamiltonians are “typical” in the sense that the statistical properties of
physical observables are the same as those predicted by the microcanonical
ensemble. [37, 38, 39]

Our quantum typicality refers to an idea that anticipates that almost all
quantum systems will have similar dynamical properties [40].

###  III.3 Random states

Controlling the dynamics of an extensive system is practically impossible when
the system becomes large and complex. Formally our control strategy requires
$N$ states Eq. (15) and (16). Can we reduce the number of states to $L<K<N$ ?
We know from Ref. [41] that sampling quantum states at random can be seen to
be induced by the sampling of unitaries at random. Sampling a set of
unitaries, which are bounded, means we sample a set of actions we can apply to
our state, which is global.

We will now sample using a random state, as follows:

$|\Psi({\vec{\theta}})\rangle=\frac{1}{\sqrt{N}}\sum_{j}^{N}e^{i\theta_{j}}\left|\phi_{j}\right\rangle$
(21)

where $\vec{\theta}=(\theta_{1},\theta_{2}...\theta_{N})$ is a vector of
random phases and $N$ is the size of the Hilbert space [42, 43]. Employing a
random set defined by different $\vec{\theta}$ we can resolve the identity:

$\hat{I}=\lim_{M\rightarrow\infty}\frac{1}{M}\sum_{l}|\psi(\vec{\theta}_{l})\rangle\langle\psi(\vec{\theta}_{l})|$
(22)

Employing now $K$ random states to sample the control and relying on quantum
typicality we expect $K<N$.

##  IV Results

Our primary goal in this study is to develop a method to control systems with
large Hilbert space while minimizing the computational effort. The objective
of the present model was to cool by increasing the system’s purity at the
final steady-state after multiple cycles using the OCT algorithm. It is
important to note that this objective is a multi-state problem. We have shown
that such a task is possible [22], but the main drawback is the computation
scaling of the problem with the system’s size.

Assuming that the target of the process is to reach a colder state, any
proposed mechanism for the process has to maintain the population of the
target state while allowing the population of all other states to repopulate
selectively. After many excitation-relaxation cycles, the memory of the
initial state will be erased, and the obtained transformation will be of which
mentioned in Eq. (1). In this model, rotational cooling of AlF Sec. II we
define the cost function $J$ to populate states which will through spontaneous
emission will populate lower $j$ values and penalize increase in $j$ values.
To control this type of system we have created several realizations for an
increasing number of random states $L$ to control the transformation. We
wanted to check:

  1. 1.

Does the total transformation can be represented.

  2. 2.

If it does, how many states $L$ are required to converge the entire system and
the cooling transformation sufficiently? (where $K$ is the size of the cooling
transformation).

The solution we obtained from the OC algorithm (the control field) was
employed as a pilot field for the following realization. We hoped to decrease
the complexity of the problem by reducing the amount of iteration needed to
solve the control problem for the next realization.

Employing the control field from a small $L$, as a pilot field, to a larger
one did not significantly improve the convergence. Such a behavior has been
observed for generic quantum objectives [44].

Figure 1 displays the control fidelity of systems with a different number of
random states $L$ as a function to the number of iterations. To estimate the
number of states needed to converge the control, we have sampled a small set
and checked the obtained control field on the full-sized (Hilbert space)
dynamics.

We have a controlled scheme of a single field polarized to the $Z$ spatial
axis. The initial population distributed among different random states
corresponds to a system with $J=11\to N=430$. The target was to move the
population to all states in Hilbert space, up to excited state at $J\leqslant
10$ and eliminate any transitions to higher $J$ states that open decay
channels to even higher states. The control field can manipulate the system
towards the target state with a fitness of $99\%$. However, this fidelity fits
only the contracted space; later on, we will show that the same field acting
on the full space system results in lower fidelity. Thus, some impurities are
leaking in.

Figure 1: (main) The controlled Infidelity (log scale), of different number
random phase sampling, as a function of the number of iterations required to
converge the control($\approx{99\%}$ at the target state). The number of RP
wavefunctions is assigned by different colors in the legend. The number of
states of each model influences the number of iteration, increase
monotonically. Inset The numerical effort (log scale) as a function of the
number of RP states. The numerical effort is a function of Hilbert space times
the number of iterations, showing a polynomial scaling.

In the inset 1 we show how the Numerical effort increases polynomially. The
effort is defined by the following: ${N_{Effort}}=\ln({N_{I}}\times{M_{RP}})$
where $M_{RP}$ is the number of random phase states used in this model and
$N_{I}$ states the number of iterations needed for the OC algorithm to
converge to the predetermined thresh hold.

As mentioned, we have used the controlled field on the full system; figure 2
shows the infidelity of the entire system as a function of the number of
random phase states $L$. To check the quality of the transition we used the
field from the controlled scheme on the dynamics of the entire system. We find
a monotonic growth of fidelity corresponding to the number $L$ of random
states.

Figure 2:  The infidelity as a function of the number of random phase base
functions $L$. The infidelity is calculated for the full transformation.

The fidelity of the full transformation as expected is smaller than the one
obtained from the control calculated from a finite set of random phase states.
To find the number of random states which lead to convergence, we extrapolate
our results, taking its infidelity as a measure.

Figure 3:  Fitting the infidelity of the full transformation (red dots) as a
function of the number of RP states (blck curve). Extrapolating to high
fidelity allows to estimate the effective size of the cooling transformation
(44,0.01).

In figure 3 we have fitted our calculated points to an exponential curve.
Defining an acceptable threshold we get the number of sampling states $K$
required to converge the full system of dimension $N$.

The marked point is a good guess to the complexity of the transformation
corresponding to an effective transformation size of $K=44$. Nevertheless, the
computation effort for this sampling space is vary high. The effective
transformation size is still much smaller than the dimension of the full
transformation.

A concrete measure for cooling in our context, is the change in normalized
entropy. Employing Eq. (14) we define the normalized entropy decrease:

$\Delta{S_{eff}}=\frac{{S_{FS}^{RP}-S_{initial}^{J=11}}}{{S_{Th}^{J=10}-S_{initial}^{J=11}}}$
(23)

where $S_{FS}^{RP}$ is the entropy of the full system obtained from the
control sequence, $S_{initial}^{J=11}$ is the thermal entropy of the initial
state (J=11), and $S_{Th}^{J=10}$ is the thermal entropy of target for our
transformation. The normalization is with respect to the thermal entropy
difference between $J=11$ to $J=10$ (at T= $\sim 30^{0}K$).

Figure 4:  The monotonic decrease of normalized entropy, Eq. (23) with $L$,
for the complete transformation obtained by each RP model (blue stars). The
entropy value of the thermal state at J=10 is marked by a red line and for J=9
by a blue line. The objective of full transformation was to cool from $J=11$
to J=10, but it is clear that for the 19-RP model we accomplished an even
colder state.

Figure 4 shows the effective change in entropy for each model (RP). A
monotonic decrease of the effective entropy is a clear indication for cooling.

The samples describe different realisation, where the initial state were
created randomly. We have used these realisations for the same transformation
and used it on the full spaced system. The fidelity we have got on the full
system with respect to the desired transformation as calculated, and built for
the calculation of the standard deviation, which is defined by
$STD=\frac{\sqrt{<J^{2}>-<J>^{2}}}{|<J>|}$. In the next figure we show how the
STD behaves in different models.

Figure 5:  The Standard Deviation (STD) of the infidelity of the total
transformation as a function of the number of samples ($n$). The blue curve is
data obtained by optimizing a single RP model. The orange curve was obtained
for simultaneous transformation for 3-RP sampling. The red curve was obtained
for simultaneous transformation for 6-RP sampling.

An ensemble of typical states should have a small standard deviation with
respect to a local observable. This is confirmed by Figure 5 displaying the
standard deviation of the target infidelity when the size of the sample
increases. Each sample is composed of independent random realizations with
different $L$. We expect convergence as $\frac{1}{\sqrt{n}}$, where $n$ is the
sample size. The standard deviation is quite small for the random phase
wavefunction sampling. As expected when random transformation $L$ is larger
($3-RP$ and $6-RP$) the STD decreases.

##  V conclusions

Laser cooling of the internal degrees of freedom of a molecule is a difficult
task due to the large occupied Hilbert space. A shaped pulse generates a
unitary transformation accompanied by spontaneous emission. We optimize the
unitary transition ${\cal U}$ such as after many cooling cycles the target
state is colder than the initial one. We have utilized quantum typicality for
computing the cooling map. The transformation is calculated with an increasing
number of random phase wavefunctions. We show convergence (of infidelity) to
the target state. Doing so while reducing the computational effort.

We have shown that this method can be utilized to model cooling of internal
degrees of freedom of molecules. It is anticipated that a series of such
transformation can cool the rotation to its ground state.

##  VI Acknowledgment

We thank Gerard Meijer for his helpful discussions. Work supported by Israel
Science Foundation.

## References

  * [1] S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, et al., “Training schrödinger’s cat: quantum optimal control,” The European Physical Journal D, vol. 69, no. 12, p. 279, 2015.
  * [2] G. M. Huang, T. J. Tarn, and J. W. Clark, “On the controllability of quantum-mechanical systems,” Journal of Mathematical Physics, vol. 24, no. 11, pp. 2608–2618, 1983.
  * [3] F. Albertini and D. D’Alessandro, “Notions of controllability for bilinear multilevel quantum systems,” IEEE Transactions on Automatic Control, vol. 48, no. 8, pp. 1399–1403, 2003.
  * [4] V. Ramakrishna, M. V. Salapaka, M. Dahleh, H. Rabitz, and A. Peirce, “Controllability of molecular systems,” Physical Review A, vol. 51, no. 2, p. 960, 1995.
  * [5] A. P. Peirce, M. A. Dahleh, and H. Rabitz, “Optimal control of quantum-mechanical systems: Existence, numerical approximation, and applications,” Physical Review A, vol. 37, no. 12, p. 4950, 1988.
  * [6] R. Kosloff, S. A. Rice, P. Gaspard, S. Tersigni, and D. Tannor, “Wavepacket dancing: Achieving chemical selectivity by shaping light pulses,” Chemical Physics, vol. 139, no. 1, pp. 201–220, 1989.
  * [7] J. P. Palao and R. Kosloff, “Optimal control theory for unitary transformations,” Physical Review A, vol. 68, no. 6, p. 062308, 2003.
  * [8] I. Serban, J. Werschnik, and E. K. U. Gross, “”Optimal control of time-dependent targets”,” vol. 71, p. 053810, 2005.
  * [9] J. Werschnik, and E. K. U. Gross, “”Quantum optimal control theory”,” J. Phys. B., vol. 40, pp. R175–R211, SEP 28 2007.
  * [10] R. Eitan, M. Mundt, and D. J. Tannor, “Optimal control with accelerated convergence: Combining the krotov and quasi-newton methods,” Physical Review A, vol. 83, no. 5, p. 053426, 2011.
  * [11] D. M. Reich, M. Ndong, and C. P. Koch, “Monotonically convergent optimization in quantum control using krotov’s method,” The Journal of chemical physics, vol. 136, no. 10, p. 104103, 2012.
  * [12] R. Kosloff, “Time-dependent quantum-mechanical methods for molecular dynamics,” The Journal of Physical Chemistry, vol. 92, no. 8, pp. 2087–2100, 1988.
  * [13] S. Arora and B. Barak, Computational complexity: a modern approach.  Cambridge University Press, 2009.
  * [14] C. Bartsch and J. Gemmer, “Dynamical typicality of quantum expectation values,” Physical review letters, vol. 102, no. 11, p. 110403, 2009.
  * [15] J. Von Neumann, Mathematical foundations of quantum mechanics.  Princeton university press, 2018.
  * [16] R. Uzdin and S. Rahav, “Global passivity in microscopic thermodynamics,” Physical Review X, vol. 8, no. 2, p. 021064, 2018.
  * [17] R. Krems, B. Friedrich, and W. C. Stwalley, Cold molecules: theory, experiment, applications.  CRC press, 2009.
  * [18] J. Carini, S. Kallush, R. Kosloff, and P. Gould, “Enhancement of ultracold molecule formation using shaped nanosecond frequency chirps,” Physical review letters, vol. 115, no. 17, p. 173003, 2015.
  * [19] S. Kallush and R. Kosloff, “Unitary photoassociation: One-step production of ground-state bound molecules,” Physical Review A, vol. 77, no. 2, p. 023421, 2008.
  * [20] M. Yeo, M. T. Hummon, A. L. Collopy, B. Yan, B. Hemmerling, E. Chae, J. M. Doyle, and J. Ye, “Rotational state microwave mixing for laser cooling of complex diatomic molecules,” Physical review letters, vol. 114, no. 22, p. 223003, 2015.
  * [21] C.-Y. Lien, C. M. Seck, Y.-W. Lin, J. H. Nguyen, D. A. Tabor, and B. C. Odom, “Broadband optical cooling of molecular rotors from room temperature to the ground state,” Nature communications, vol. 5, no. 1, pp. 1–7, 2014.
  * [22] A. Aroch, S. Kallush, and R. Kosloff, “Optimizing the multicycle subrotational internal cooling of diatomic molecules,” Physical Review A, vol. 97, no. 5, p. 053405, 2018.
  * [23] W. M. Itano, J. C. Bergquist, J. J. Bollinger, and D. J. Wineland, “Cooling methods in ion traps,” Physica Scripta, vol. 1995, no. T59, p. 106, 1995\.
  * [24] M. Zeppenfeld, B. G. Englert, R. Glöckner, A. Prehn, M. Mielenz, C. Sommer, L. D. van Buuren, M. Motsch, and G. Rempe, “Sisyphus cooling of electrically trapped polyatomic molecules,” Nature, vol. 491, no. 7425, pp. 570–573, 2012.
  * [25] Y. Maday and G. Turinici, “New formulations of monotonically convergent quantum control algorithms,” The Journal of Chemical Physics, vol. 118, no. 18, pp. 8191–8196, 2003.
  * [26] D. M. Reich and C. P. Koch, “Cooling molecular vibrations with shaped laser pulses: optimal control theory exploiting the timescale separation between coherent excitation and spontaneous emission,” New Journal of Physics, vol. 15, no. 12, p. 125028, 2013.
  * [27] T. Brixner, G. Krampert, T. Pfeifer, R. Selle, G. Gerber, M. Wollenhaupt, O. Graefe, C. Horn, D. Liese, and T. Baumert, “Quantum control by ultrafast polarization shaping,” Physical review letters, vol. 92, no. 20, p. 208301, 2004.
  * [28] S. Truppe, S. Marx, S. Kray, M. Doppelbauer, S. Hofsaess, H. C. Schewe, B. Sartakov, and G. Meijer, “Spectroscopy, buffer gas cooling and radiation pressure slowing of alf molecules,” Bulletin of the American Physical Society, vol. 64, 2019.
  * [29] S. Truppe, S. Marx, S. Kray, M. Doppelbauer, S. Hofsäss, H. C. Schewe, N. Walter, J. Pérez-Ríos, B. G. Sartakov, and G. Meijer, “Spectroscopic characterization of aluminum monofluoride with relevance to laser cooling and trapping,” Physical Review A, vol. 100, no. 5, p. 052513, 2019.
  * [30] R. N. Zare, Angular momentum, vol. 33.  Wiley, New York, 1988.
  * [31] A. Lenard, “Thermodynamical proof of the gibbs formula for elementary quantum systems,” Journal of Statistical Physics, vol. 19, no. 6, pp. 575–586, 1978.
  * [32] D. J. Tannor, Introduction to quantum mechanics: a time-dependent perspective.  2007\.
  * [33] J. Werschnik and E. Gross, “Quantum optimal control theory,” Journal of Physics B: Atomic, Molecular and Optical Physics, vol. 40, no. 18, p. R175, 2007\.
  * [34] I. Degani, A. Zanna, L. SÆlen, and R. Nepstad, “Quantum control with piecewise constant control functions,” SIAM Journal on Scientific Computing, vol. 31, no. 5, pp. 3566–3594, 2009.
  * [35] J. P. Palao and R. Kosloff, “Quantum computing by an optimal control algorithm for unitary transformations,” Physical review letters, vol. 89, no. 18, p. 188301, 2002.
  * [36] D. J. Tannor, V. Kazakov, and V. Orlov, “Control of photochemical branching: Novel procedures for finding optimal pulses and global upper bounds,” in Time-dependent quantum molecular dynamics, pp. 347–360, Springer, 1992\.
  * [37] M. Rigol, V. Dunjko, and M. Olshanii, “Thermalization and its mechanism for generic isolated quantum systems,” Nature, vol. 452, no. 7189, pp. 854–858, 2008.
  * [38] J. von Neumann, “Proof of the ergodic theorem and the h-theorem in quantum mechanics,” The European Physical Journal H, vol. 35, no. 2, pp. 201–237, 2010.
  * [39] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, “From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics,” Advances in Physics, vol. 65, no. 3, pp. 239–362, 2016.
  * [40] P. Reimann, “Typicality for generalized microcanonical ensembles,” Physical review letters, vol. 99, no. 16, p. 160404, 2007.
  * [41] P. F. Romero, Equilibration and Typicality in Quantum Processes.  PhD thesis, Monash University, 2020.
  * [42] M. Nest and R. Kosloff, “Quantum dynamical treatment of inelastic scattering of atoms at a surface at finite temperature: The random phase thermal wave function approach,” The Journal of chemical physics, vol. 127, no. 13, p. 134711, 2007.
  * [43] S. Kallush and S. Fleischer, “Orientation dynamics of asymmetric rotors using random phase wave functions,” Physical Review A, vol. 91, no. 6, p. 063420, 2015.
  * [44] S. Kallush and R. Kosloff, “Scaling the robustness of the solutions for quantum controllable problems,” Physical Review A, vol. 83, no. 6, p. 063412, 2011.


# Short-term forecasting COVID-19 cumulative confirmed cases: Perspectives for
Brazil

Matheus Henrique Dal Molin Ribeiro<EMAIL_ADDRESS>Ramon
Gomes da Silva  Viviana Cocco Mariani  Leandro dos Santos Coelho  Industrial &
Systems Engineering Graduate Program (PPGEPS), Pontifical Catholic University
of Parana (PUCPR). 1155, Rua Imaculada Conceicao, Curitiba, PR, Brazil.
80215-901  Department of Mathematics, Federal Technological University of
Parana (UTFPR). Via do Conhecimento, KM 01 - Fraron, Pato Branco, PR, Brazil.
85503–390  Mechanical Engineering Graduate Program (PPGEM), Pontifical
Catholic University of Parana (PUCPR). 1155, Rua Imaculada Conceicao,
Curitiba, PR, Brazil. 80215-901  Department of Electrical Engineering, Federal
University of Parana (UFPR). 100, Avenida Coronel Francisco Heraclito dos
Santos, Curitiba, PR, Brazil. 81530-000

###### Abstract

The new Coronavirus (COVID-19) is an emerging disease responsible for
infecting millions of people since the first notification until nowadays.
Developing efficient short-term forecasting models allow knowing the number of
future cases. In this context, it is possible to develop strategic planning in
the public health system to avoid deaths. In this paper, autoregressive
integrated moving average (ARIMA), cubist (CUBIST), random forest (RF), ridge
regression (RIDGE), support vector regression (SVR), and stacking-ensemble
learning are evaluated in the task of time series forecasting with one, three,
and six-days ahead the COVID-19 cumulative confirmed cases in ten Brazilian
states with a high daily incidence. In the stacking learning approach, the
cubist, RF, RIDGE, and SVR models are adopted as base-learners and Gaussian
process (GP) as meta-learner. The models’ effectiveness is evaluated based on
the improvement index, mean absolute error, and symmetric mean absolute
percentage error criteria. In most of the cases, the SVR and stacking ensemble
learning reach a better performance regarding adopted criteria than compared
models. In general, the developed models can generate accurate forecasting,
achieving errors in a range of 0.87% - 3.51%, 1.02% - 5.63%, and 0.95% - 6.90%
in one, three, and six-days-ahead, respectively. The ranking of models in all
scenarios is SVR, stacking ensemble learning, ARIMA, CUBIST, RIDGE, and RF
models. The use of evaluated models is recommended to forecasting and monitor
the ongoing growth of COVID-19 cases, once these models can assist the
managers in the decision-making support systems.

###### keywords:

ARIMA , COVID-19 , Forecasting , Decision-making , Machine learning , Time-
series

††journal: Chaos, Solitons & Fractals

##  1 Introduction

The new Coronavirus (COVID-19) is an emerging disease responsible for
infecting millions of people and killing thousands worldwide since the first
notification until nowadays, according to the World Health Organization (WHO)
[1, 2]. Also according to WHO, Brazil registered 40.581 confirmed cases until
April 22th 2020, holding the 12th position in the world ranking in the number
of confirmed cases of COVID-19, and 2nd position in the Americas (behind the
United States of America).

Due to the impacts of the COVID-19 pandemic in people’s lives and the world’s
economy, the governments and population are most concerned with (i) when the
COVID-19 outbreak will peak; (ii) how long the outbreak will last and (iii)
how many people will eventually be infected [3]. Further, Boccaletti et al.
[4] have identified at least three scientific communities that may cooperate
in the effort to deal with the current pandemic: (i) the community of applied
mathematicians, virologists and epidemiologists, developing sophisticated
diffusion models to the specific properties of a given pathogen; (ii) the
community of complex systems scientists who study the spread of infections
using compartmental models, using methods and principles from statistical
mechanics and nonlinear dynamics; and (iii) the community of scientists who
incorporate artificial intelligence (AI) and most specifically deep learning
approaches to produce accurate predictive models. Also, different studies are
evaluating the impacts of COVID-19 on society, whether through predictions of
future cases, as well as variables capable of helping to understand the spread
of this disease [5, 6, 7, 8, 9].

Moreover, epidemiological time series forecasting plays an important role in
health public system, once it allows the managers to develop strategic
planning to avoid possible epidemics. Forecasting diseases as accurate as
possible is important due to their impact on the public health system. To
ensure this accuracy, AI models have been widely used to forecast
epidemiological time series over the years [10, 11, 12]. Moreover, in the AI
context, Vaishya et al. [13] presented a review of trends in COVID-19 data
analysis.

Regarding this context, the objective of this paper is to explore and compare
the predictive capacity of machine learning regression and statistical models,
in the task of forecasting one, three, and six-days-ahead COVID-19 cumulative
cases in Brazil. In this respect, datasets of ten Brazilian states some with a
high incidence of COVID-19 until now, like Sao Paulo and Rio de Janeiro, are
adopted to evaluates the forecasting efficiency through of the autoregressive
integrated moving average (ARIMA), cubist regression (CUBIST), random forest
(RF), ridge regression (RIDGE), support vector regression (SVR), and stacking-
ensemble learning models. In the stacking learning modelling, which is an
effective ensemble learning approach [14, 15], CUBIST, RF, RIDGE, and SVR are
used as base-learners (weak models), and Gaussian process (GP) as meta-learner
(strong model). The out-of-sample forecasting accuracy of each model is
compared by some performance metrics such as the improvement percentage index
(IP), mean absolute errors (MAE), and symmetric mean absolute percentage error
(sMAPE).

The contributions of this paper can be summarized as follows:

  * •

The first contribution is related to the presentation of a novel analysis of
the forecast model for cumulative confirmed cases of COVID-19 in Brazil, whose
accuracy of the models assists governors in decision-making to contain the
pandemic and strategies concerning the health system;

  * •

The second contribution, we can highlight the use of heterogeneous machine
learning models, as well as the stacking-ensemble learning approach to
forecast the Brazilian cumulative confirmed cases of COVID-19;

  * •

Also, this paper evaluates models forecasting in a multi-day-ahead forecasting
strategy. The forecasting time horizons are the interval of one, three, and
six-days-ahead. This range of the forecasting time horizon allows us to verify
the effectiveness of the predicting models in different scenarios, helping in
future strategies in fighting COVID-19.

The remainder of this paper is organized as follows: Section 2.1 a brief
description of the dataset adopted in this paper is given. The forecasting
models applied in this study are described in Section 2.2. Section 3 details
the procedures applied in the research methodology. Results obtained and
related discussion about models’ forecasting performance are given on Section
4. Finally, Section 5 concludes this work with considerations and some
directions for future research proposals.

##  2 Material and Methods

This section presents the description of the material analyzed (Section 2.1)
as well as the models description applied in this paper (Section 2.2).

###  2.1 Dataset Description

The collected dataset refers to the cumulative confirmed cases of COVID-19
that occurred in Brazil until April, 18 or 19 of 2020. The dataset was
collected from an API [16] that retrieves the daily information about COVID-19
cases from all 27 Brazilian State Health Offices, gather them, and make it a
publicly available. Among the 27 federative units (26 states and one federal
district), ten states some with a high incidence of COVID-19 cases and other
states with lower temperatures, states from south of Brazil, were chosen,
among them are Amazonas (AM), Bahia (BA), Ceara (CE), Minas Gerais (MG),
Parana (PR), Rio de Janeiro (RJ), Rio Grande do Norte (RN), Rio Grande do Sul
(RS), Santa Catarina (SC), and Sao Paulo (SP). The measurement period of each
state varies, once each state counts since the day of its first case until the
day of the last report. The cumulative confirmed cases and deaths of each
state, as well as the period from the first and last reports, are illustrated
in Table 1. The change in the way of accounting for the number of cases, by
the health departments, may change the data presented here.

Table 1: First and last report dates by state State |  | Number of
---
observed days
First report | Last report |  | Cumulative
---
confirmed cases
| Cumulative
---
deaths
AM | 34 | 13/03/2020 | 19/04/2020 | 2044 | 182
BA | 43 | 06/03/2020 | 19/04/2020 | 1249 | 45
CE | 35 | 16/03/2020 | 19/04/2020 | 3306 | 189
MG | 42 | 08/03/2020 | 19/04/2020 | 1154 | 39
PR | 36 | 12/03/2020 | 18/04/2020 | 960 | 49
RJ | 38 | 05/03/2020 | 19/04/2020 | 4675 | 402
RN | 30 | 12/03/2020 | 18/04/2020 | 561 | 26
RS | 38 | 10/03/2020 | 19/04/2020 | 869 | 26
SC | 39 | 12/03/2020 | 19/04/2020 | 1025 | 35
SP | 53 | 25/02/2020 | 19/04/2020 | 14267 | 1015

A heatmap of the cumulative confirmed cases is presented in Figure 1.

Figure 1: Heatmap of the cumulative confirmed cases of the analyzed states

###  2.2 Methodologies

This section describes a brief of each model employed in the data analysis.

  * •

ARIMA is a Box & Jenkins modelling usually employed to deal with non-
stationary time series. In fact, the ARIMA model is full specified by
autoregressive (p), different degrees of trend differences (d), and moving
average operators (q). These parameters are used do define the model order,
and usually defined by grid-search, as well as by autocorrelation and partial
autocorrelation function. In this context, the model is described as
ARIMA(p,d,q) [17].

  * •

CUBIST is a rule-based model, which performs predictions following the
regression of trees principle [18]. Through the use of a committee of the
rules, and using the neighborhood concept similar to $k$-nearest-neighbor
modelling, the final forecasting is obtained.

  * •

GP is composed of a set of random variables Gaussian distributed and fully
specified by its mean and covariance (kernel) function [19]. In this paper,
the GP with a linear kernel is adopted.

  * •

RIDGE is a regularized regression approach [20] which employs a penalization
term in the ordinary least squares algorithm. It is an effective tool, once it
reduces the bias of parameter estimates by controlling the standard errors.
Moreover, the model can deal with inputs multi-collinearity problem.

  * •

RF is a bagging ensemble-based model, which combines the bagging advantages
characterized by the creation of multiple samples, with refitting through of
the bootstrap technique, from the same set of data, and random selection of
predictors to compose each node of the decision tree [21]. RF is a fast and
robust supervised learning method able to deal with the randomness of the time
series. Furthermore, it is interesting because, in addition to being an
ensemble approach, only the number of predictors for each node needs to be
tuned.

  * •

SVR consists in determining support vectors (points) close to a hyperplane
that maximizes the margin between two-point classes obtained from the
difference between the target value and a threshold. To deal with non-linear
problems SVR takes into account kernel functions, which calculates the
similarity between two observations. In this paper, the linear kernel is
adopted. The main advantages of the use of SVR lies in its capacity to capture
the predictor non-linearity and then use it to improve the forecasting cases.
In the same direction, it is advantageous to employ this perspective in this
case study adopted, since that the samples are small [22].

  * •

Stacked Generalization or stacking is an ensemble-based approach [23] which
combines through a meta-learner the predictions of a set of weak models (base-
learners) to obtain a stronger learner. This approach usually operates into
two levels, where in the first level the base-learners are trained and its
predictions are obtained. In the next stage, a meta-learner uses, as inputs,
the predictions of the previous level in the training phase. The stacking
predictions are obtained from meta-learner. The main advantage of the stacking
ensemble is that this approach can improve the accuracy and additionally
reduce error variance [14].

##  3 Proposed forecasting framework

This section describes the main steps in the data analysis adopted by CUBIST,
RF, RIDGE, SVR, and stacking models. Also, the ARIMA modelling is described.

Step 1: Firstly, the raw data is split into training and test datasets. The
test dataset is composed of six last observations, and the training dataset by
the remain samples [14]. The training data are centered by its mean value and
divided by its standard deviation. To develops multi-days-ahead COVID-19 cases
forecasting, recursive strategy is employed [24]. In this aspect, one model is
fitted for one-day-ahead forecasting. Next, the recursive strategy uses the
forecasting value as an input for the same model to forecast the next step,
continuing this manner until reaching the desirable horizon. The training
structure adopted in this paper is stated as follows,

$\displaystyle
y_{(t+1)}=f\left\\{y_{t},\ldots,y_{t+1-n_{y}}\right\\}+\epsilon\quad\epsilon\sim
N(\mathbf{0},\sigma^{2}),$ (1)

in which $f$ is a function related to the adopted model in the training stage,
$y_{t+1}$ is the COVID-19 case one-day-ahead, $n_{y}=5$ are the past confirmed
cases, $\epsilon$ is the random error, following a normal distribution with
zero mean and constant variance. In this paper, the aim is to obtain the cases
up to $H$ next days, especially up to 1 (ODA, one-day-ahead), 3 (TDA, three-
days-ahead), and 6-days-ahead (SDA, six-days-ahead), respectively. The
following structures are considered,

$\hat{y}_{t+h}=\begin{cases}f\left[y_{t},y_{t-1},\ldots,y_{t},\ldots,y_{t-n_{y}+1}\right]&\text{if
}h=1\\\
f\left[\hat{y}_{t+h-1},\ldots,\hat{y}_{t+h},\ldots,y_{t},\ldots,y_{t+h-n_{y}}\right]&\text{if
}h\in[2,n_{y}]\\\
f\left[\hat{y}_{t+h-1},\ldots,\hat{y}_{t+h-n_{y}}\right]&\text{if
}h\in[n_{y}+1,\ldots,H],\end{cases}$ (2)

where $\hat{y}_{t+h}$ is the forecast value at time $t$ and forecast horizon
up to $h$, $y_{t+h-n_{y}}$ and $\hat{y}_{t+h-n_{y}}$ are the previously
observed and forecast cases lags in $n_{y}=5$ days. The $n_{y}$ value is
chosen through grid-search with purpose to capture the best data behavior.

Step 2: In the stacking modelling, the base-learners CUBIST, RF, RIDGE, SVR
are trained and its forecasting are used as inputs for meta-learner GP. In the
training stage, leave-one-out cross-validation with a time slice is adopted
[14]. Finally, the out-of-sample forecasts are computed. These approaches are
developed using the caret package [25]. The ARIMA modeling is performed
through the use of forecast package [26, 27] with use of auto.arima function.
To define the ARIMA order, grid-search is adopted, and the most suitable order
is that reach a lower Akaike and Bayesian Akaike criteria information. Both
analyses are developed using R software [28]. All hyperparameters employed in
this study are presented in Table 1 in B.

Step 3: To evaluate the effectiveness of adopted models, from obtained
forecasts out-of-sample (test set), performance IP (3), MAE (4), and sMAPE (5)
criteria are computed as

$\displaystyle\operatorname{MAE}$ $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left|y_{i}-\hat{y}_{i}\right|,$ (3)

$\displaystyle\operatorname{sMAPE}$ $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left|\frac{\hat{y}_{i}-y_{i}}{(|y_{i}|+|\hat{y}_{i}|/2)}\right|,$
(4) $\displaystyle\operatorname{IP}$ $\displaystyle=$ $\displaystyle
100\times\frac{\displaystyle M_{c}-M_{b}}{M_{c}},$ (5)

where $n$ is the number of observations, $y_{i}$ and $\hat{y}_{i}$ are the
i-th observed and predicted values, respectively. Also, the $M_{c}$ and
$M_{b}$ represent the performance measure of compared and best models,
respectively.

Figure 2 presents the proposed forecasting framework.

Figure 2: Proposed forecasting framework

##  4 Results

This section describes the results of the developed experiments in forecasts
out-of-sample (test set). First, Section 4.1 compares the results of evaluated
models over ten datasets and three forecasting horizons adopted. In Table 1 in
A, the best results regarding accuracy are presented in bold. Additionally,
Figures 3 up to 4 illustrate the relation between observed and predicted
values achieved by models with best set of performance measures depicted in
Table 1, as well as box-plots for out-of-sample errors are illustrated in
Figure 5.

###  4.1 Performance Measures for compared models

In this section, the main results achieved by the best model regarding MAE and
sMAPE criteria are presented for short-term forecasting multi-days- ahead of
cumulative cases of COVID-19 from ten Brazilian states.

  * •

AM: In this state, CUBIST, and RIDGE approaches could be considered to
forecasting COVID-19 cases. In fact, in respect to ODA and TDA, CUBIST
outperforms models, while for SDA the RIDGE achieves better accuracy regarding
MAE and sMAPE than others. The improvement in the MAE for ODA and TDA achieved
by CUBIST ranges between 6.58% - 92.77%, and 11.39% - 88.54%, respectively.
Through sMAPE analysis, the RIDGE model outperforms other models, and this
criterion is reduced in the range of 16.46% - 91.88%, for SDA horizon.

  * •

BA, MG, RS, and SP: For these states, in all forecasting windows, the SVR
approach achieved better accuracy than other models, for both MAE and sMAPE
criteria in the multi-days-ahead forecasting task of the confirmed number of
COVID-19. In fact, the improvement in sMAPE is ranged in 13.26% - 95.11%,
4.23% - 94.88%, and 38.59% - 95.24%, respectively, in ODA, TDA, and SDA
forecasting horizons. Moreover, the same behavior is observed when the
improvement in sMAPE criterion is obtained.

  * •

CE and RN: In the CE state, the ARIMA model has a better performance in the
forecasting out-of-sample than other models for ODA and TDA time windows. In
this aspect, for MAE criterion, the improvement is ranged between 72.36% -
98.03%, and 45.93% - 92.40%, for ODA, and TDA time windows, respectively. For
sMAPE, the improvement on ODA, and TDA horizons is 65.06% - 97.84%, and 32.81%
- 92.53%, respectively. The SVR has better results than ARIMA model for SDA.
Considering the RN state, the same analysis is developed for ODA, and TDA
horizons. The exception to the SDA horizon, in which the CUBIST model has
better effectiveness in the MAE and sMAPE criteria than remain models.

  * •

PR, RJ, and SC: For these states localized into the south region (PR and SC)
and southeast region (RJ) of Brazil, the most appropriate approach to forecast
cumulative cases of COVID-19 is the stacking ensemble, exception in ODA
horizon, when ARIMA model has better results. Stacking overcomes the drawback
of single models and achieves the best accuracy than other models. In fact,
for these states, the improvement in MAE and sMAPE are between 14.01% -
94.68%, and 17.48% - 95.41%, respectively, for ODA horizon. The improvement in
order forecasting horizons presents the same behavior of ODA, with the
greatest magnitude of improvement for TDA and SDA.

Remark: In this experiment, 180 scenarios (10 datasets, 3 forecasting
horizons, and 6 models) were evaluated for the task of forecasting cumulative
COVID-19 cases. In an overview, the best models for each state, obtained sMAPE
ranged between 0.87% - 3.51%, 1.02% - 5.63%, and 0.95% - 6.90% for ODA, TDA,
and SDA forecasting, respectively. The ranking of models in all scenarios is
SVR, stacking ensemble, ARIMA, CUBIST, RIDGE, and RF models. In contrast to
finds of [29], for the datasets evaluated in this paper, ARIMA modelling was
effective in some situation for very-short horizons When the horizon is SDA,
ARIMA model has worst performance than most of compared models. However, for
ODA the applications are limited. From a broader perspective, the efficiency
of SVR is due to its ability to deal with small size dataset, while the
stacking ensemble combines the advantages of several single models to learn
the data behavior and obtain forecasts similar to observed values. On the
other hand, the difficulty of the RF model to forecasting cumulative COVID-19
cases could be attributed to the fact that this approach requires more
observations to effectively learn the data pattern.

(a) AM  (b) BA

(c) CE  (d) MG

Figure 3: Predicted versus observed cumulative confirmed cases of COVID-19 for
AM, BA, CE, and MG states

(a) PR  (b) RJ

(c) RN  (d) RS

(e) SC  (f) SP

Figure 4: Predicted versus observed cumulative confirmed cases of COVID-19 for
PR, RJ, RN, RS, SC, and SP states

According to the information depicted in Figures 3 and 4 it is possible to
identify that the behavior of the data is learned by the evaluated models,
which can forecasting compatible cases with the observed values. The good
performance obtained in the training phase persists in the test stage. In the
Figures 3a and 4c the models, RIDGE and CUBIST, as well as in Figures 3d and
4f SVR presented difficulties to capture the variability of the first
observations. The dataset is reduced for all states, which justifies the
difficulties of the mathematical models to learn the behavior.

Figure 5 shows the box-plots of out-of-sample forecasting errors in the SDA
horizon for each model and dataset used. This horizon is chosen to analysis
due to the recursive strategy adopted, once the errors increase according to
the growth of the forecasting horizon. The box diagram depicts the variation
of absolute errors for each model, which reflects the stability of each model.
In this context, the dots out of boxes are considered outliers errors, and the
black dot inside of the box is the MAE for each model.

Figure 5: Box-plot for absolute error according to model and state for
COVID-19 forecasting up to SDA

Through the box-plot analysis, boxes with lower size indicate models with
lower variation in the errors, and the results presented in Table 1 are
corroborated by the depicted in Figure 5. Models with lower errors also reach
better stability, which means that the most suitable modelling for each state
can maintain a learning pattern, achieving homogeneous prediction errors.

##  5 Conclusion and Future Research

In this paper, six machine learning approaches named CUBIST, RF, RIDGE, SVR,
and stacking ensemble, as well as ARIMA statistical model, were employed in
the task of forecasting one, three, and six-days-ahead the COVID-19 cumulative
confirmed cases in ten Brazilian states with a high daily incidence. The
COVID-19 cumulative confirmed cases for AM, BA, CE, MG, PR, RJ, RN, RS, SC,
and SP states were used. The IP, MAE, and sMAPE criteria were adopted to
evaluate the performance of the compared approaches. Moreover, the stability
of out-of-sample errors was evaluated through box-plots.

In respect of obtained results, it is possible to infer that SVR and stacking-
ensemble learning model are suitable tools to forecast COVID-19 cases for most
of the adopted states, once that these approaches were able to learn the
nonlinearities inherent to the evaluated epidemiological time series. Also,
ARIMA can be considered in some aspects for ODA, while CUBIST and RIDGE models
deserve attention for the development of this task in TDA and SDA time
windows. Therefore, the ranking of models in all scenarios is SVR, stacking
ensemble, ARIMA, CUBIST, RIDGE, and RF models. However, even though the models
discussed in this paper presented forecasting cases similar to those observed,
they should be used cautiously. This fact is attributed to the chaotic
dynamics of the analyzed data, as well as the diversity of exogenous factors
that can affect the daily notifications of COVID-19.

For future works, it is intended (i) to adopt deep learning approaches
combined to stacking ensemble, (ii) to employ copulas functions for data
augmentation dealing with small samples, (iii) to use multi-objective
optimization to tune hyperparameters of adopted models, (iv) to adopt set of
features which can help to explain the future cases of the COVID-19.

## CRediT Author Statement

Matheus Henrique Dal Molin Ribeiro: Conceptualization, Methodology, Formal
analysis, Validation, Writing - Original Draft, Writing - Review & Editing.
Ramon Gomes da Silva: Conceptualization, Methodology, Formal analysis,
Validation, Writing - Original Draft, Writing - Review & Editing. Viviana
Cocco Mariani: Conceptualization, Writing - Review & Editing. Leandro dos
Santos Coelho: Conceptualization, Writing - Review & Editing.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgments

The authors would like to thank the National Council of Scientific and
Technologic Development of Brazil – CNPq (Grants number: 307958/2019-1-PQ,
307966/2019-4-PQ, 404659/2016-0-Univ, 405101/2016-3-Univ), PRONEX ‘Fundação
Araucária’ 042/2018, and Coordenação de Aperfeiçoamento de Pessoal de Nível
Superior - Brasil (CAPES) - Finance Code 001 for financial support of this
work.

## References

  * World Health Organization (2020) [WHO] World Health Organization (WHO), .  Coronavirus(COVID-19).  2020\.  URL: https://covid19.who.int/; (accessed in 22 April, 2020).
  * Sohrabi et al. [2020] Sohrabi, C., Alsafi, Z., O’Neill, N., Khan, M., Kerwan, A., Al-Jabir, A., et al.  World health organization declares global emergency: A review of the 2019 novel coronavirus (COVID-19).  International Journal of Surgery 2020;76:71–76.  doi:10.1016/j.ijsu.2020.02.034.
  * Zhang et al. [2020] Zhang, X., Ma, R., Wang, L..  Predicting turning point, duration and attack rate of COVID-19 outbreaks in major western countries.  Chaos, Solitons & Fractals 2020;:109829doi:10.1016/j.chaos.2020.109829.
  * Boccaletti et al. [2020] Boccaletti, S., Ditto, W., Mindlin, G., Atangana, A..  Modeling and forecasting of epidemic spreading: The case of Covid-19 and beyond.  Chaos, Solitons & Fractals 2020;doi:10.1016/j.chaos.2020.109794.
  * Becerra et al. [2020] Becerra, M., Jerez, A., Aballay, B., Garcés, H.O., Fuentes, A..  Forecasting emergency admissions due to respiratory diseases in high variability scenarios using time series: A case study in chile.  Science of The Total Environment 2020;706:134978.  doi:10.1016/j.scitotenv.2019.134978.
  * Fanelli and Piazza [2020] Fanelli, D., Piazza, F..  Analysis and forecast of COVID-19 spreading in china, italy and france.  Chaos, Solitons & Fractals 2020;134:109761.  doi:10.1016/j.chaos.2020.109761.
  * Fong et al. [2020] Fong, S.J., Li, G., Dey, N., Crespo, R.G., Herrera-Viedma, E..  Composite monte carlo decision making under high uncertainty of novel coronavirus epidemic using hybridized deep learning and fuzzy rule induction.  Applied Soft Computing 2020;doi:10.1016/j.asoc.2020.106282.
  * Roosa et al. [2020] Roosa, K., Lee, Y., Luo, R., Kirpich, A., Rothenberg, R., Hyman, J., et al.  Real-time forecasts of the COVID-19 epidemic in China from february 5th to february 24th, 2020.  Infectious Disease Modelling 2020;5:256 – 263.  doi:10.1016/j.idm.2020.02.002.
  * Effenberger et al. [2020] Effenberger, M., Kronbichler, A., Shin, J.I., Mayer, G., Tilg, H., Perco, P..  Association of the COVID-19 pandemic with internet search volumes: A Google trendstm analysis.  International Journal of Infectious Diseases 2020;doi:10.1016/j.ijid.2020.04.033.
  * Davis et al. [2019] Davis, J.K., Gebrehiwot, T., Worku, M., Awoke, W., Mihretie, A., Nekorchuk, D., et al.  A genetic algorithm for identifying spatially-varying environmental drivers in a malaria time series model.  Environmental Modelling & Software 2019;119:275–284.  doi:10.1016/j.envsoft.2019.06.010.
  * Ribeiro et al. [2019] Ribeiro, M.H.D.M., da Silva, R.G., Fraccanabbia, N., Mariani, V.C., Coelho, L.d.S..  Forecasting epidemiological time series based on decomposition and optimization approaches.  In: 14th Brazilian Computational Intelligence Meeting (CBIC). Belém, Brazil; 2019, p. 1–8.
  * Scavuzzo et al. [2018] Scavuzzo, J.M., Trucco, F., Espinosa, M., Tauro, C.B., Abril, M., Scavuzzo, C.M., et al.  Modeling dengue vector population using remotely sensed data and machine learning.  Acta Tropica 2018;185:167 – 175.  doi:10.1016/j.actatropica.2018.05.003.
  * Vaishya et al. [2020] Vaishya, R., Javaid, M., Khan, I.H., Haleem, A..  Artificial intelligence (ai) applications for covid-19 pandemic.  Diabetes & Metabolic Syndrome: Clinical Research & Reviews 2020;14(4):337–339.  doi:10.1016/j.dsx.2020.04.012.
  * Ribeiro and Coelho [2020] Ribeiro, M.H.D.M., Coelho, L.d.S..  Ensemble approach based on bagging, boosting and stacking for short-term prediction in agribusiness time series.  Applied Soft Computing 2020;86(105837).  doi:10.1016/j.asoc.2019.105837.
  * Moreno et al. [2019] Moreno, S.R., da Silva, R.G., Ribeiro, M.H.D.M., Fraccanabbia, N., Mariani, V.C., Coelho, L.d.S..  Very short-term wind energy forecasting based on stacking ensemble.  In: 14th Brazilian Computational Intelligence Meeting (CBIC). Belém, Brazil; 2019, p. 1–8.
  * Justen [2020] Justen, A..  Covid-19: Coronavirus newsletters and cases by municipality per day.  2020\.  URL: https://brasil.io/api/dataset/covid19/caso/data/?place_type=state; (accessed in 20 April, 2020).
  * Box et al. [2015] Box, G.E., Jenkins, G.M., Reinsel, G.C., Ljung, G.M..  Time series analysis: forecasting and control.  5 ed.; John Wiley & Sons; 2015\.  ISBN 978-1-118-67502-1.
  * Quinlan [1993] Quinlan, J.R..  Combining instance-based and model-based learning.  In: Proceedings of the Tenth International Conference on International Conference on Machine Learning. ICML’93; San Francisco, CA, USA: Morgan Kaufmann Publishers Inc.  ISBN 1-55860-307-7; 1993, p. 236–243.
  * Rasmussen [2004] Rasmussen, C.E..  Gaussian Processes in Machine Learning; chap. 4.  Heidelberg, Germany: Springer.  ISBN 978-3-540-28650-9; 2004, p. 63–71.
  * Hoerl and Kennard [1970] Hoerl, A.E., Kennard, R.W..  Ridge regression: Biased estimation for nonorthogonal problems.  Technometrics 1970;12(1):55–67.  doi:10.1080/00401706.1970.10488634.
  * Breiman [2001] Breiman, L..  Random forests.  Machine Learning 2001;45(1):5–32.  doi:10.1023/A:1010933404324.
  * Drucker et al. [1997] Drucker, H., Burges, C.J.C., Kaufman, L., Smola, A.J., Vapnik, V..  Support vector regression machines.  In: Mozer, M.C., Jordan, M.I., Petsche, T., editors. Advances in Neural Information Processing Systems 9. MIT Press; 1997, p. 155–161.
  * Wolpert [1992] Wolpert, D.H..  Stacked generalization.  Neural Networks 1992;5(2):241–259.  doi:10.1016/S0893-6080(05)80023-1.
  * Moreno et al. [2020] Moreno, S.R., da Silva, R.G., Mariani, V.C., Coelho, L.d.S..  Multi-step wind speed forecasting based on hybrid multi-stage decomposition model and long short-term memory neural network.  Energy Conversion and Management 2020;213(112869).  doi:10.1016/j.enconman.2020.112869.
  * Kuhn [2008] Kuhn, M..  Building predictive models in R using the Caret package.  Journal of Statistical Software, Articles 2008;28(5):1–26.  doi:10.18637/jss.v028.i05.
  * Hyndman et al. [2020] Hyndman, R., Athanasopoulos, G., Bergmeir, C., Caceres, G., Chhay, L., O’Hara-Wild, M., et al.  forecast: Forecasting functions for time series and linear models; 2020.  URL: http://pkg.robjhyndman.com/forecast; R package version 8.12.
  * Hyndman and Khandakar [2008] Hyndman, R.J., Khandakar, Y..  Automatic time series forecasting: the forecast package for R.  Journal of Statistical Software 2008;26(3):1–22.  URL: http://www.jstatsoft.org/article/view/v027i03.
  * R Core Team [2018] R Core Team, .  R: A language and environment for statistical computing.  R Foundation for Statistical Computing; Vienna, Austria; 2018.
  * Benvenuto et al. [2020] Benvenuto, D., Giovanetti, M., Vassallo, L., Angeletti, S., Ciccozzi, M..  Application of the arima model on the covid-2019 epidemic dataset.  Data in Brief 2020;29(105340).  doi:10.1016/j.dib.2020.105340.

##  Appendix A Performance Measures

Table 1 presents the performance measures for each model in each state and
forecasting horizon.

Table 1: Performance measures for each evaluated model

State |  Forecasting Horizon | Criteria | Model
---|---|---|---
ARIMA | CUBIST | RF | RIDGE | Stacking | SVR
AM | ODA | MAE | 95 | 45 | 622.17 | 48.17 | 121.5 | 56.33
sMAPE | 6.61% | 2.80% | 42.50% | 2.83% | 7.13% | 3.18%
TDA | MAE | 101.33 | 71.33 | 622.17 | 83.67 | 176.67 | 80.5
sMAPE | 6.55% | 4.50% | 42.50% | 4.49% | 10.47% | 4.19%
SDA | MAE | 119.17 | 162.17 | 622.17 | 62.33 | 233.17 | 79.17
sMAPE | 6.97% | 9.55% | 42.50% | 3.45% | 13.87% | 4.13%
BA | ODA | MAE | 12 | 93.83 | 366.33 | 45.33 | 107.67 | 42.33
sMAPE | 1.56% | 9.16% | 42.02% | 4.36% | 10.68% | 4.15%
TDA | MAE | 70 | 132 | 366.33 | 74.33 | 171.67 | 59.67
sMAPE | 8.00% | 12.92% | 42.02% | 7.46% | 17.32% | 5.63%
SDA | MAE | 155.67 | 152.33 | 366.33 | 152.83 | 215.83 | 73.17
sMAPE | 15.41% | 15.08% | 42.02% | 15.16% | 22.25% | 6.90%
CE | ODA | MAE | 18 | 65.17 | 916 | 70.33 | 220.83 | 87.67
sMAPE | 0.87% | 2.49% | 40.28% | 2.81% | 8.20% | 3.17%
TDA | MAE | 69.66 | 128.83 | 916 | 149.83 | 382.17 | 136.67
sMAPE | 3.01% | 4.48% | 40.28% | 5.39% | 14.48% | 4.78%
SDA | MAE | 257 | 118.17 | 916 | 98.17 | 484.33 | 164.17
sMAPE | 9.34% | 4.11% | 40.28% | 3.52% | 18.78% | 5.77%
MG | ODA | MAE | 32 | 17.5 | 235.5 | 24.33 | 56.5 | 16
sMAPE | 3.63% | 1.81% | 26.21% | 2.50% | 5.59% | 1.57%
TDA | MAE | 26 | 21.33 | 235.5 | 21.67 | 78.17 | 21
sMAPE | 3.08% | 2.20% | 26.21% | 2.13% | 7.81% | 2.04%
SDA | MAE | 55 | 36.83 | 235.5 | 32.17 | 97.83 | 14.33
sMAPE | 5.43% | 3.58% | 26.21% | 3.14% | 9.88% | 1.41%
PR | ODA | MAE | 31 | 27.33 | 163.5 | 38 | 23.5 | 35.33
sMAPE | 3.96% | 3.26% | 21.09% | 4.50% | 2.69% | 4.18%
TDA | MAE | 51.66 | 57.33 | 163.5 | 76.5 | 28.17 | 60.17
sMAPE | 6.21% | 6.56% | 21.09% | 8.61% | 3.21% | 6.89%
SDA | MAE | 73.67 | 118 | 163.5 | 151 | 24.17 | 117.17
sMAPE | 8.20% | 12.56% | 21.09% | 15.75% | 2.75% | 12.53%
RJ | ODA | MAE | 110 | 165.5 | 1305.67 | 273.67 | 69.5 | 360.83
sMAPE | 3.17% | 3.82% | 37.06% | 6.25% | 1.70% | 8.09%
TDA | MAE | 120 | 275.67 | 1305.67 | 462.83 | 68 | 429.33
sMAPE | 3.18% | 6.24% | 37.06% | 10.20% | 1.65% | 9.49%
SDA | MAE | 158.33 | 532.67 | 1305.67 | 696.17 | 65.17 | 529.5
sMAPE | 3.67% | 11.34% | 37.06% | 14.67% | 1.58% | 11.43%
RN | ODA | MAE | 6 | 17 | 152.5 | 24.83 | 30.33 | 18.33
sMAPE | 1.61% | 3.87% | 39.28% | 5.56% | 6.45% | 4.14%
TDA | MAE | 8.33 | 30.83 | 152.5 | 37.67 | 54 | 35.5
sMAPE | 2.11% | 6.54% | 39.28% | 8.51% | 11.66% | 7.69%
SDA | MAE | 36.33 | 15.83 | 152.5 | 62 | 54 | 18.5
sMAPE | 7.61% |  3.42%% | 39.28% | 12.76% | 11.66% | 4.15%
RS | ODA | MAE | 12 | 12.83 | 146.67 | 11.33 | 45.5 | 8.17
sMAPE | 1.64% | 1.62% | 19.82% | 1.43% | 5.76% | 0.97%
TDA | MAE | 24 | 19.17 | 147.33 | 18.67 | 71.33 | 8.5
sMAPE | 3.22% | 2.47% | 19.92% | 2.42% | 9.14% | 1.02%
SDA | MAE | 34.5 | 34.17 | 147.5 | 37.67 | 91.83 | 7.83
sMAPE | 4.31% | 4.26% | 19.95% | 4.74% | 11.89% | 0.95%
SC | ODA | MAE | 21 | 93.67 | 179.5 | 180.5 | 33.83 | 177.67
sMAPE | 2.43% | 9.66% | 20.97% | 17.53% | 3.66% | 17.27%
TDA | MAE | 44.33 | 100.33 | 179.5 | 277 | 41 | 257.33
sMAPE | 4.76% | 10.30% | 20.97% | 25.34% | 4.39% | 23.79%
SDA | MAE | 56 | 102.83 | 179.5 | 338.5 | 43.83 | 330.33
sMAPE | 5.65% | 10.53% | 20.97% | 29.95% | 4.68% | 29.23%
SP | ODA | MAE | 436 | 1587 | 3799 | 537.33 | 1363.83 | 409
sMAPE | 4.65% | 13.47% | 35.85% | 4.44% | 11.44% | 3.51%
TDA | MAE | 1485.66 | 2471.83 | 3801 | 579.17 | 2243 | 326.67
sMAPE | 14.56% | 21.81% | 35.88% | 4.79% | 19.47% | 2.77%
SDA | MAE | 2779 | 3054.67 | 3801.5 | 591.83 | 2665.83 | 362.83
sMAPE | 24.74% | 27.60% | 35.88% | 4.95% | 23.55% | 3.04%

##  Appendix B Hyperparameters

Table 1 presents the hyperparameters obtained by grid-search for the models
employed in this paper. In the stacking modeling, for the GP meta-learner
there is no hyperaparameter to be tuned.

Table 1: Hyperparameters selected by grid-search for each evaluated model

State | Model
---|---
ARIMA | CUBIST | SVR | RIDGE | RF
(p,d,q) | Committees | Neighbors | Cost | Regularization |  | Number of randomly
---
selected predictors
AM | (1,2,0) | 10 | 5 | 1 | 3.16E-03 | 2
BA | (0,2,1) | 20 | 9 | 1 | 1E-04 | 2
CE | (2,2,1) | 1 | 9 | 1 | 0 | 4
MG | (0,2,1) | 1 | 9 | 1 | 1E-04 | 2
PR | (0,2,1) | 20 | 5 | 1 | 3.16E-03 | 3
RJ | (0,2,1) | 1 | 9 | 1 | 1E-04 | 3
RN | (1,1,0) | 1 | 9 | 1 | 3.16E-03 | 5
RS | (0,1,0) | 1 | 9 | 1 | 1E-04 | 3
SC | (0,2,1) | 10 | 0 | 1 | 3.16E-03 | 5
SP | (0,2,0) | 20 | 9 | 1 | 1E-04 | 5