<figure><img src="image_1.png" /></figure>

Article

# Data-Driven Framework for Uncovering Hidden Control Strategies in Evolutionary Analysis

Nourddine Azzaoui <fnref n="1" />,*,†, Tomoko Matsui <fnref n="2" />,†, and Daisuke Murakami <fnref n="2" />

**Abstract:** We devised a data-driven framework for uncovering hidden control strategies used by an evolutionary system described by an evolutionary probability distribution. This innovative framework enables deciphering of the concealed mechanisms that contribute to the progression or mitigation of such situations as the spread of COVID-19. Novel algorithms are used to estimate the optimal control in tandem with the parameters for evolution in general dynamical systems, thereby extending the concept of model predictive control. This marks a significant departure from conventional control methods, which require knowledge of the system to manipulate its evolution and of the controller's strategy or parameters. We use a generalized additive model, supplemented by extensive statistical testing, to identify a set of predictor covariates closely linked to the control. Using real-world COVID-19 data, we delineate the descriptive behaviors of the COVID-19 epidemics in five prefectures in Japan and nine countries. We compare these nine countries and group them on the basis of shared profiles, providing valuable insights into their pandemic responses. Our findings underscore the potential of our framework as a powerful tool for understanding and managing complex evolutionary processes.

**Keywords:** data-driven optimization algorithm; model predictive control; evolutionary probability distribution; generalized additive model; classification; COVID-19 evolution

## 1. Introduction

The evolution of natural dynamical systems is a complex process where the dynamics are influenced by a multitude of factors, including epidemic evolution, environmental changes, and various exogenous random events. The traditional methods of studying evolutionary processes often involve creating simplified models that may not fully capture this complexity. When it comes to influencing or controlling this evolution, optimal control theory (OCT) is a powerful and flexible tool. This provides valuable insights into the likely paths of evolution and the factors that may affect it. The aim of this work is to investigate an evolutionary process, such as epidemic evolution, from the perspective of OCT in order to uncover the hidden mechanisms that lead to deterioration or improvement of the situation. We assume that the evolutionary process is described by a controlled system via a hidden mechanism to which we do not have access and that we cannot change directly; we can only observe its effects through measurements related to its evolution. We also assume that this hidden control mechanism is affected by an infinite number of random parameters related to the behavior and characteristics of societies and populations. We start from a simple linear control model based on Kalman-type filtering. We aim to jointly determine the parameters of the spatial-temporal process and the active optimal control that leads to the observed evolution. The determination of the appropriate parameters is based on the best proximity between the observed data and that generated by the simple model.

<figure><img src="image_2.png" /></figure>

### check for updates

**Citation:** Azzaoui, N.; Matsui, T.; Murakami, D. Data-Driven Framework for Uncovering Hidden Control Strategies in Evolutionary Analysis. *Math. Comput. Appl.* **2023**, *28*, 103. https://doi.org/10.3390/mca28050103

Academic Editor: Cristiana João
Soares da Silva

Received: 22 August 2023

Revised: 4 October 2023

Accepted: 13 October 2023

Published: 20 October 2023

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Laboratoire de Mathématiques Blaise Pascal, University of Clermont Auvergne, CNRS, 63000 Clermont-Ferrand, France</footnote>

<footnote n="2">² Department of Statistical Modeling, Institute of Statistical Mathematics, 10-3 Midori-cho, Tachikawa 190-8562, Tokyo, Japan; tmatsui@ism.ac.jp (T.M.); dmuraka@ism.ac.jp (D.M.)</footnote>

<footnote>* Correspondence: nourddine.azzaoui@uca.fr; Tel.: +33-473407081</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
Given that the objects in evolution are probability distributions, various customized
distance measures are employed to evaluate their proximity. The Hellinger distance is,
thus, preferred due to its superior numerical performance, inherited from its link to the
ℓ₂-norm. The coronavirus disease 2019 (COVID-19) is being examined using the principles
of OCT. For instance, Péni et al. suggest an optimal control strategy based on model
predictive control (MPC) theory for a compartmental epidemic model that depicts the
dynamics of COVID-19 virus transmission. Péni and Szederkényi explore the application
of the MPC approach to a nonlinear compartmental model. The COVID-19 outbreak
in Germany is being studied using the epidemiological SIDARTH model to determine
the best feedback control approach using MPC. Carli et al. present a unique feedback
control mechanism using an epidemiological SIR (susceptible, infectious, and recovered or
removed)-based model and the MPC approach to help policymakers effectively mitigate
the effects of the COVID-19 pandemic in numerous regions, including Italy, where the
pandemic reaches its peak. The SIR model is being explored by Morato et al. in relation
to MPC strategies to control COVID-19 infections in Brazil. All of these studies use MPC
to simulate epidemiological models such as the SIR model. However, using simulation
to investigate the COVID-19 pandemic is challenging due to the simplistic assumptions
and limitations of the models, such as their inability to represent treatments, their fixed
parameters, and their uncertainty resulting from data quality concerns. We, thus, aim to
identify the variables affecting the spread of COVID-19 by using a data-driven model and
a simple controlled Kalman-type dynamical system.

Our goal in this work is to make a simple assumption about the dynamics of COVID-
19 virus transmission and then utilize the data to determine the parameters governing
it. An extensive initial analysis of Japanese data is being performed to understand the
dynamics of the epidemic in diverse prefectures (Tokyo, Osaka, Hokkaido, Fukuoka,
and Okinawa). Their differing characteristics are compared and linked via a generalized
additive model (GAM) to different socio-demographic structures and the related policies
that are adopted in each prefecture.

After validating our data-driven model on the Japanese data, we examine whether
it can be applied to other countries or geographical regions with different population
structures and dynamics. We do this by applying an MPC-like procedure, as we do
for the Japan data, to simultaneously obtain the dynamics and control parameters that
contribute to the observed evolution of the epidemic in each country. We then perform a
comparative analysis between the countries based on their management strategies, such as
for vaccinations and lockdowns.

* Descriptive behavior of epidemic by country: Using the epidemic's fitted controlling parameters for each country or region, we attempt to clarify the relationship between the parameters and various socio-demographic indicators, as well as the various policies that governments implemented to reduce the epidemic's spread.

* Comparative analysis and clustering: Using the contrasting characteristics and indicators of the various countries, we group them on the basis of the similarities in their profiles. We then look at the key variables that led to countries being grouped together. We also conduct additional investigations for each group to determine which factors contributed most to the grouping.

Our contributions are summarized below:

* We propose using a hidden controller framework in which the optimal control and the parameters for evolution in general dynamical systems, including the COVID-19 epidemic dynamical system, are simultaneously estimated. This transforms the linear evolution scheme into a nonlinear problem.

* We present an efficient optimization algorithm from the practical point of view that is based on a classical optimization algorithm, e.g., fmincon in MATLAB.
* We present the results of a comparative analysis between prefectures and regions of Japan with differing population structures and epidemic management policies that were obtained using our hidden controller framework.

* We show that generalizing the model using our hidden controller framework to other countries (Australia, Germany, Italy, Brazil, South Africa, and so on) can reveal the effect of different policies and population structures on the evolution of the epidemic.

We discuss in detail the methodology we used, as well as the related work in Section 2. In Section 3, we describe our experiments, along with the application of our proposed model to real-world datasets. In Section 4, we discuss our results and the importance of our study. We conclude in Section 5 with a summary of the key points.

## 2. Methods

The most commonly used data-driven control technique is based on MPC (Figure 1). It enables the controller or user to interact with a dynamical system and affect its evolution by applying a specific control variable at selected times or epochs depending on the predicted horizon. Our approach was inspired by this technique but differs in two key aspects: we do not have direct access to the system, so we cannot manipulate its evolution, and we do not have any insight into the controller's strategy or parameters. (Here we focus on a Kalman-type dynamical system, but our approach can be applied to any parametric/non-parametric dynamical system). Our knowledge is limited to the fact that the controller or the system uses a model from a specific system family (here, for example, one in the form of Equation (2)). That is, we have access to only the observed final output data. This idea is sketched in Figure 2.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** In the MPC framework, users can access and manipulate (i.e., change or control) the system parameters.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** In our framework, users do not have access to the items in the shaded area; they know only the shapes of the optimizer, model, and system and are unable to act on that information.</figcaption></figure>

This section is organized as follows; We begin by reviewing the related work and classical optimal control techniques in Sections 2.1 and 2.2. In Sections 2.3 and 2.6, we describe our proposed hidden control framework.
## 2.1. Targeted Control Problem and Its Formalism

We assume a discrete-time context in which decision-makers (e.g., government agencies) observe the distribution of the disease variable (e.g., the numbers of weekly infected, recovered, and deceased individuals), possibly with observational inaccuracy, for a given time period (day, week, month, etc.). For each sub-population $P_i$, the observed distribution is generally estimated empirically from the data, which we denote $\hat{\Pi}_i(k) = (\hat{\pi}_i^1(k), \hat{\pi}_i^2(k), \dots, \hat{\pi}_i^{n+1}(k))$. The estimated evolutionary distributions for the whole population can be represented as

$$ \mathcal{M}_k = \begin{pmatrix} \hat{\pi}_1^1(k) & \dots & \hat{\pi}_1^1(k) \\ \hat{\pi}_1^2(k) & \dots & \hat{\pi}_1^2(k) \\ \vdots & \ddots & \vdots \\ \hat{\pi}_1^{n+1}(k) & \dots & \hat{\pi}_1^{n+1}(k) \end{pmatrix}. \quad (1) $$

Since $\pi_i^{n+1}(k)$ can be expressed as $\pi_i^{n+1}(k) = 1 - \sum_{l=1}^{n} \pi_i^l(k)$, we are interested in the dynamical evolution of the discrete $n$-dimensional probability vector $x_i(k) = (\pi_i^1(k), \pi_i^2(k), \dots, \pi_i^n(k))$. Because we estimate the optimal control by sub-population, we omit the suffix $i$ hereafter to simplify the description (e.g., $x_i(k)$ is represented as $x_k$ and a column of $\mathcal{M}_k$ as $M_k$). This simplified description represents the state vector of the disease distribution at time $k$. To model the evolution, we propose using a discrete Kalman-like controlled system:

$$ x_{k+1} = A x_k + B u_k, \quad \forall k = 1, \dots, K-1, \quad (2) $$

where

* The term $u_k \in \mathbb{R}^m$ is an $m$-dimensional vector describing the unknown optimal control. This unknown control is the hidden mechanism that led to the observed shape of the evolution.

* The deterministic state matrix $A = (a_{ij})_{i,j=1,...,n}$, defined in $\mathbb{R}^{n \times n}$ and describing the logical evolution of $x_k$ when no control is applied, is also unknown.

* The unknown control matrix $B = (b_{ij})_{i=1,...,n, j=1,...,m}$, defined in $\mathbb{R}^{n \times m}$, describes how the control is applied to the system to bring it to its targeted state.

The idea is to use only observed data $\mathcal{M}_k$, in Equation (1), to simultaneously estimate all the unknown objects of this control problem, namely, matrices $A, B$, and the control term $u_k$, which leads to the observed shape of the data. When $A$ and $B$ are known, the optimal control $\hat{u}_k^*$ can be obtained using classical techniques that have been widely discussed in the literature (e.g., [1–3]). In our case, these parameters are unknown and should be estimated from the observed data. Our strategy is to use the optimal control results obtained using an iterative algorithm that minimizes a cost function.

## 2.2. Control Algorithm and Related Work

For given deterministic matrices $A$ and $B$, the solution of the discrete-time optimal control problem, Equation (2), of a linear system is described in [1,3]. It consists of using a linear optimization scheme, Equation (2), and finding the optimal control $u = (u_1, \dots, u_K)$ that minimizes the objective or cost function:

$$ J(u) = x'_{K} S x_{K} + \sum_{k=1}^{K-1} x'_{k} Q x_{k} + u'_{k} R u_{k}, \quad (3) $$

where exponent $(\cdot)'$ represents the transpose, and $Q, R$, and $S$ are cost matrices defined as follows.

* $Q \in \mathbb{R}^{n \times n}$ is a positive definite weighting matrix associated with the state vector. It is linked to the relative importance assigned to state vector $x_k$ in cost function $J(\cdot)$.
In other words, a higher value in the diagonal of this matrix indicates a bigger penalty or greater emphasis placed on minimizing the related state vector $x_k$.

* $R \in \mathbb{R}^{m \times m}$ is a positive definite weighting matrix affecting the relative relevance or weight assigned to control $u_k$ in cost function $J(\cdot)$. It specifies how much emphasis the control algorithm should place on minimizing the control effort. A higher value in its diagonal, similar to that in the Q matrix, indicates a higher penalty or greater importance placed on reducing the control input.

* $S$ is a positive definite weighing matrix linked to the final state vector. It is a cross-weighting matrix used by cost function $J(\cdot)$ to penalize discrepancies between the inputs from state $x_k$ and control $u_k$. It is aimed at identifying any cross-dependencies between the states and control inputs. We are not interested in the final state here, so we suppose that $S = 0$.

The optimal control problem has been widely studied (e.g., [1,2]). In our case, a linear quadratic regulator variant is used to minimize, Equation (3). The solution of the state space, Equation (2), of a discrete time-invariant linear system is given in a more general framework [4] for $2 \le k \le K$:

$$x_k = A^k x_1 + \sum_{i=1}^{K} A^{k-i} B u_i. \quad (4)$$

This equation expresses state vector $x_k$ in terms of the initial state vector $x_1$ and a sequence of control vectors, $u_i, i = 1, \dots, k$, that are implicitly dependent: $A, B$, and weighting matrices $Q, R$. In this setting, the controlled system, Equation (2), with dimension $K$ is controllable if the controllability matrix,

$$C = [B, AB, \dots, A^i B, \dots, A^{K-1} B], \quad (5)$$

has column rank $K$ (i.e., $K$ linearly independent columns). Under these conditions on parameters $A, B$, the optimal control solution [4] is given by

$$u_k^* = -K x_k. \quad (6)$$

In this setting, the closed-form expression of the gain matrix $K$ is obtained using a linear quadratic regulator (LQR). This is performed iteratively in two steps using a backward approach, as detailed in Algorithm 1:

* Gain matrix $K$ is calculated from evolutionary parameter matrices $A, B, Q$, and $R$ via $S_{k+1}$ at step $k$:

$$K(A, B) = (R + B'S_{k+1}B)^{-1}B'S_{k+1}A. \quad (7)$$

* $S_k$ is obtained from the solution of the Riccati equation [5]:

$$S_k = Q + A'S_{k+1}A - A'S_{k+1}BK(A, B). \quad (8)$$

By combining Equation (2) and (6), we can write the evolutionary structure of the controlled system in Equation (2) recursively:

$$x_k^* = (A - B K(A, B)) x_{k-1}. \quad (9)$$

In the general framework of optimal control, the parameters are known and fixed. In our study, we are dealing with unknown parameters. Therefore, our objective becomes the estimation of these unidentified matrices in the Kalman filter parameters, namely, $A, B$, as well as weighting matrices $Q, R$.
**Algorithm 1:** Computing optimal control and its corresponding trajectory (see Figure 3)

**Require:** Estimated distributions $M_k$ in Equation (1) from observed data.

**Require:** A vector $\mathcal{Y}$ of dimension $n(n+m+2)$ of the system parameters.

1. Retrieve $A, B, Q, R$ from $\mathcal{Y}$.

**Ensure:** Matrices $A, B$ verify controllability condition (5).

**Ensure:** $Q, R$ are symmetric positive definite.

2. Initialize $S_{K+1}$.

3. **for** $k = K - 1$ to 0 **do**

4. Compute gain $\mathcal{K}(A, B) \leftarrow (R + B'S_{K+1}B)^{-1}B'S_{K+1}A$ as in Equation (7).

5. Update $S_k \leftarrow Q + A'S_{K+1}A - A'S_{K+1}B\mathcal{K}(A, B)$.

6. **end for**

7. **for** $k = 0$ to $K - 1$ **do**

8. Get input $u_k^* = -\mathcal{K}(A, B)x_k$.

9. Update using (2) to get new estimate $x_{k+1}^*$.

10. **end for**

11. **Output:** Gain matrix $\mathcal{K}(A, B)$.

12. **Output:** $x^* = (x_1^*, \dots, x_K^*)$, the controlled trajectory.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Structure of main algorithm for estimating hidden control process.</figcaption></figure>

## 2.3. Our Hidden Controller Framework

Instead of the traditional objects found in control problems, our methodology relies on using discrete probability distributions as the dynamic system object, as shown in Equation (2), to characterize the system’s behavior. A further distinguishing characteristic of our methodology is its ability to estimate both the parameters of the model and their corresponding optimal control, relying exclusively on observed data. To determine unknown parameters $A$ and $B$, our approach utilizes two interrelated optimization steps:

1. The first step is to find optimal control features $\mathcal{K}(A, B)$ and $x_k^*$, as shown in Equation (9) under controllability conditions on $A$ and $B$.

2. The second step is to minimize a suitable distance that was adapted to probability distributions between $x_k^*$ and observation $x_k$ under two constraints: $x_k$ is a discrete probability distribution calculated from observed data $M_k$, and the system is controllable. This leads to the following minimization problem:

$$ (\hat{A}, \hat{B}) = \underset{A, B, \text{s.t. } \forall t, 0 < \sum_{k=1}^{K} (A-B\mathcal{K})x_{t-1} < 1}{\operatorname{argmin}} \sum_{k=1}^{K} d(x_k^*, x_k). \quad (10) $$

The constrained minimum is taken over all $A$ and $B$ such that the system is controllable, as detailed in Algorithm 2.
**Algorithm 2:** Estimating hidden control parameters and trajectory

**Require:** Estimated distributions $M_k$ in (1) from observed data.

1. Define a convergence threshold $\epsilon$.

2. Initialize $A, B, Q, R$.

3. Transform $Y = (a_{11}, \dots, a_{nn}, b_{11}, \dots, b_{nm}, \alpha_1, \dots, \alpha_n, \beta_1, \dots, \beta_n)$.

4. Define objective function $F(.)$ in (12) for given distance.

**Ensure:** Constraint (C1) and (C2) on $Y$.

5. **while** ($F(Y) > \epsilon$) **do**

6. Call Algorithm 1 with $Y$, and get $\mathcal{K}(A, B)$ and $x^* = (x_1^*, \dots, x_K^*)$.

7. Solve quadratic programming subproblem under constraints $C_1$ and $C_2$.

8. Update $Y$ and calculate its corresponding $F(Y)$.

9. **end while**

10. **Output:** Optimal parameters $(\hat{A}, \hat{B}, \hat{Q}, \hat{R})$.

11. **Output:** Corresponding optimal control $\hat{u}^*$.

In accordance with the concept of the probability distribution, distance $d(.,.)$ in Equation (10) determines the degree of similarity between two probability distributions. A variety of metrics and distances that have been used in different circumstances and scenarios have been reported. We next discuss five of them.

## 2.4. Metrics and Distances

* The total variation distance measures the maximum difference between the probabilities that two distributions assign to an event. Therefore, it gives a worst-case measure of the difference between the two distributions, which makes it a conservative discrepancy measure. It is given by

$$d_{TV}(D_1, D_2) \triangleq \max_{S \subseteq D} (D_1(S) - D_2(S)) \\ = \frac{1}{2} \sum_{x \in D} |D_1(x) - D_2(x)| = \frac{1}{2} \|D_1 - D_2\|_1. \qquad (11)$$

For the discrete case, the right-hand side of Equation (11) can be expressed as the $l_1$ distance. Despite its frequent usage, the TV distance metric has certain disadvantages. It is sensitive to absolute differences between probabilities, which can lead to misleading results when applied to large event sets. Though it satisfies several metric properties such as non-negativity, symmetry, and triangle inequality, it is not a “proper” metric as it can assign a distance of zero to non-identical distributions. The TV distance also lacks the ability to capture dependencies among variables in multivariate distributions and does not provide information about the ‘location’ and ‘spread’ of differences. Furthermore, it is not a differentiable function, posing challenges for optimization algorithms. Finally, its computational complexity increases substantially for continuous distributions, especially in high-dimensional scenarios.

* Kullback-Leibler (KL) divergence quantifies the information loss when distribution $D_1$ is used to approximate $D_2$:

$$d_{KL}(D_1, D_2) = \sum_{x \in D} D_1(x) \log \left( \frac{D_1(x)}{D_2(x)} \right).$$

The KL divergence has several advantages and disadvantages when used as a measure of the difference between two probability distributions. Among its advantages are its high sensitivity to discrepancies in probability distributions, including those occurring in distribution tails, and its widespread application across multiple fields, such as in-
formation theory and machine learning. It also has a deep connection with maximum likelihood estimation, which enhances its value for theoretical analysis. The disadvantages include its lack of symmetry; that is, the divergence from distribution $D_1$ to $D_2$ is not the same as that from $D_2$ to $D_1$, which can make it counterintuitive as a “distance” measure. Furthermore, it is not a metric as it does not satisfy properties such as symmetry and the triangle inequality. Another critical disadvantage is its undefined nature for non-overlapping distributions. Finally, when probability distributions are estimated from the data, KL divergence can be excessively sensitive to the sample size, leading to potentially unreliable estimates for smaller sample sizes.

* The Wasserstein distance, also known as the earth mover's distance, measures the minimum amount of work needed to transform one pile into another pile. For continuous or discrete probability distributions $D_1$ and $D_2$, the $p$-Wasserstein distance ($p \ge 1$) is defined as

$$d_{W_p}(D_1, D_2) \triangleq \inf \left\{ \left( \mathbb{E}[|X - Y|^p] \right)^{\frac{1}{p}}; \mathbb{P}_X = D_1 \text{ and } \mathbb{P}_Y = D_2 \right\} \\ = \inf \left\{ \left( \sum |x_i - y_i|^p P(X = x_i, Y = y_i) \right)^{\frac{1}{p}} \right\}.$$

The inf bound is taken over all random variable couples $(X, Y)$ such that the distribution of $X$ is $D_1$ and the distribution of $Y$ is $D_2$. It offers an intuitive geometric interpretation, acting as a measure of the minimum cost required to transform one distribution into another. In addition, it captures shifts in location, a property not captured by other metrics, making it advantageous in cases in which location is important. Moreover, unlike measures such as KL divergence, the Wasserstein distance remains defined even for non-overlapping distributions and exhibits robustness to noise. However, its computation can be resource-intensive due to the underlying optimization problem. It, thus, poses challenges, particularly in higher dimensions. The distance also depends on the choice of a cost function, which may not always have an evident best choice and can affect the distance calculation. Although sensitivity to location shifts can be beneficial, in instances where location is irrelevant, this sensitivity could capture unnecessary differences. Finally, as many other measures, the Wasserstein metric faces scalability issues with high-dimensional data, demanding dimension reduction techniques to be effective.

* The Hellinger distance measures the total difference between the square roots of the probabilities in the two distributions:

$$d_H(D_1, D_2) = \frac{1}{2} \sum_{x \in D} \left( \sqrt{D_1(x)} - \sqrt{D_2(x)} \right)^2.$$

The Hellinger distance has several advantages. It possesses metric properties including symmetry and the satisfaction of the triangle inequality, making it a true metric. It is less sensitive to outliers due to the comparison of square roots of probabilities instead of the probabilities themselves. Its fixed range, typically between 0 and 1, aids in comparing distances across various distribution pairs. Moreover, its simplicity when comparing Gaussian distributions offers computational efficiency. However, the Hellinger distance is not without drawbacks. Its structure, rooted in square roots of probabilities, makes it less sensitive to differences in distribution tail behavior. It measures global dissimilarity, often failing to capture local differences effectively. Finally, unlike KL divergence, it does not quantify information loss when one distribution approximates another, a property that may be desirable in certain contexts.

* The $\chi^2$ distance is the sum of the squared differences between the two distributions, divided by the first distribution. It is not symmetric.
$$d_{\chi^2}(D_1, D_2) = \sum_{x \in D} \frac{(D_1(x) - D_2(x))^2}{D_1(x)}$$

Each of these metrics comes with its own strengths and limitations. For instance, KL divergence ($d_{KL}$) and the chi-squared distance ($d_{\chi^2}$) are non-symmetric, making them less intuitive for certain applications. On the other hand, metrics such as the TV distance ($d_{TV}$) tend to be more conservative. In this study, our primary objective was to ensure robust convergence and computational efficiency in the optimization process. To identify the most suitable metric, we conducted comprehensive experiments testing these distances on our dataset. We encountered challenges regarding convergence and stability during the optimal control and estimation processes, as shown in Equation (10). Our experiments revealed that the Hellinger distance ($d_H$) exhibits superior performance in terms of convergence, stability, and efficiency. This can be attributed to its connection with the Euclidean $\ell_2$ norm, which provides certain computational advantages. Therefore, it is the primary distance measure we used in this study.

## 2.5. Our Parameter Estimation

To prepare for the calculation of matrix parameters A, B of the hidden closed-loop Equation (9) on the basis of observed data, we vectorize them: $(a_{11}, \dots, a_{nn}, b_{11}, \dots, b_{nm})' \in \mathbb{R}^{n^2+nm}$. Furthermore, we assume that the hidden controller uses weighting matrices Q and R to calculate the gain matrix K in Equation (7). There is no systematic way to make this choice as it depends only on the user's preferences in terms of control over parameters A and B. Most often, the weighting matrices Q and R are chosen to be diagonal. Without a loss of generality, and due to the small number of observations, we use $Q = diag(\alpha_1, \dots, \alpha_n)$ and $R = diag(\beta_1, \dots, \beta_n)$. We then rearrange all the free parameters in a vector: $Y = (a_{11}, \dots, a_{nn}, b_{11}, \dots, b_{nm}, \alpha_1, \dots, \alpha_n, \beta_1, \dots, \beta_n)' \in \mathbb{R}^{n(n+m+2)}$. Let us now consider the objective cost function:

$$F: Y \mapsto F(Y) = F(A, B, Q, R) = \sum_{k=1}^{K} d(x_k^*, M_k), \quad (12)$$

where $d(\cdot, \cdot)$ is one of the distances $d_{TV}$, $d_H$, $d_{\chi^2}$, $d_{W_p}$, or $d_{KL}$, and

* $M_k = (\hat{\pi}_i^1(k), \hat{\pi}_i^2(k), \dots, \hat{\pi}_i^n(k))' \in \mathbb{R}^n$ is the statistically estimated distribution for population segment i at observation time k. (Index i is omitted for $M_k$ and $x_k$ to simplify the description, as mentioned in Section 2.1).

* $x_k^*$ is the closed-loop state of the controlled system. The computation of vector $x_k^*$ for $k = 1, \dots, K$ is shown in Equation (9). It is implicitly dependent on the parameters in vector Y. This step is nested inside the general optimization problem, Equation (10). Vector Y should satisfy two constraints, $C_1$ and $C_2$, on the one hand, to comply with both the controllability constraints on A and B and the positive definiteness of Q and R and, on the other hand, to satisfy the probability structure of state vector $x_k$.

**C1:** Vector Y should verify the constraint that A and B lead to a controllable system by verifying Equation (5). Moreover, the unknown weighting matrices Q and R must be positive definite. For instance, in simulations and experiments, we selected Q and R as diagonal matrices, which necessitates that their diagonal elements be positive. This implies that all elements $\alpha_1, \dots, \alpha_n$ and $\beta_1, \dots, \beta_n$ are positive.

**C2:** For all $k$, vector $x_k^*$ calculated using Equation (9) should satisfy a probability structure, meaning that all its elements should be positive and their sum is less than 1.
The final optimal parameter vector, $\hat{Y}$, is then obtained by solving the nonlinear optimization problem under constraints:

$$ (\hat{A}, \hat{B}, \hat{Q}, \hat{R}) = \underset{Y \in C_1 \cap C_2}{\mathbf{argmin}} [F(Y)]. \quad (13) $$

From the practical point of view and to reduce the relative error in estimating the unknown coefficients of the parameters, we assume that the hidden controller resets vector $x_k$ to the observed value each time when calculating it for interval $K$. To implement this assumption, we divide the range of time-ordered observations $x_k$ by a pre-selected value $T$ (e.g., for the COVID-19 dataset, $T = 3$ (three days) or $T = 7$ (one week)) and denote the number of observation groups formed by $L = [\frac{K}{T}]$, where $[x]$ is the integer part. We then consider the partition group $\{G_1, \dots, G_L\}$ on all the observations $\{M_1, \dots, M_K\}$ ordered in accordance with time: the observations in group $G_1$ being $G_1 = \{M_1, \dots, M_T\}$, those in group $G_2$ being $G_2 = \{M_{T+1}, \dots, M_{2*T}\}$, and those in the last group being $G_{nstep} = \{M_{(L-1)*T+1}, \dots, M_{L*T}\}$, for example. For any $1 \le \ell \le L$, we take vector $x_k$ as follows:

* $x_k = M_k$ for $k = (\ell - 1) * T + 1$. Thus, we reset vector $x_k$ to the first observed value in each group.

* $x_k = (A + BK) * x_{k-1}$ for all $k = (\ell - 1) * T + 2, \dots, T * \ell$.

For the solution of the nonlinear optimization problem (10), we use sequential quadratic programming, which is an efficient and widely used technique [6–9], Moreover, it is readily available in the MATLAB software used for the optimization part of this work.

The results produced by Algorithm 2 are dependent on the outcomes from Algorithm 1. To secure reliability in the results, Algorithm 1 is configured to perform repeated calculations, utilizing different initial values, until it successfully converges to an optimal solution. In cases where convergence is not readily achieved, a systematic approach is employed to modify initial values, thus, facilitating the necessary conditions for convergence, and thereby ensuring the robustness of the results obtained.

These algorithms were applied to each segment of the population. The trajectories estimated using these algorithms and the observed data for different countries are plotted in Figure 5.

## 2.6. Interpretation of Control Parameters

In accordance with the preceding discussions, we computed the controlled part, denoted here as $\hat{B}, \hat{u}^*$, of the resulting dynamical system. This computation enabled us to estimate the hidden parameters and, in turn, uncover the controller’s strategy. Another key accomplishment of this study was correlating these controlling elements with the characteristics of the corresponding population. The next step is to, a posteriori, explore how these elements affected the evolution of the dynamical system. To this end, we fit the derived controlled part to the observed population characteristics as illustrated in Figure 3. For each population segment $P_i$, this was performed using a multivariate GAM. Such models are well-known for their adaptability in dealing with various forms of data and their capacity to model nonlinear relationships. We conducted extensive statistical testing on the GAM output to identify the set of variables closely linked with the control part.

### Multivariate GAM

Let us consider a multivariate GAM in which the response is vector $X = \hat{B}\hat{u}^* \in \mathbb{R}^n$ and denote $C^{(1)}, \dots, C^{(J)}$ as the dependent covariables describing the characteristics of the population. As shown in Table 1, we consider six characteristic COVID-19 variables specific to Japanese prefectures. The responses and covariates (predictors) are observed at times $k = 1 \dots K$.

$$ X = \mu_0 + S_1(C^{(1)}) + S_2(C^{(2)}) + \dots + S_J(C^{(J)}) + \varepsilon \quad (14) $$
The functions $S_1(\cdot), \dots, S_J(\cdot)$ are flexible smooth functions and can fit complex patterns in the data that would not be captured by a linear model. At the same time, they maintain interpretability: the fitted functions can be plotted and easily interpreted. The specifics of fitting a GAM depend on the type of response variable (whether it is binary, count data, continuous, etc.), the chosen link function, and the chosen form for the nonlinear functions $S_1(\cdot), \dots, S_J(\cdot)$. This model has been well studied, and many software packages are available for fitting GAMs, such as Mixed GAM Computation Vehicle with Automatic Smoothness Estimation (mgcv) in R [10,11]. As depicted in Figure 4, the application of the GAM yields an array of newly minted comprehensible features designated as $z^{(1)}, \dots, z^{(J)}$. These features serve to quantify the intensity of the association between the control and corresponding population characteristics, thereby offering a method to gauge the extent of a particular population characteristic's effect on the control of the hidden dynamical system.

<table><caption>**Table 1.** Candidate weekly predictor covariates for five Japanese prefectures.</caption><thead><tr><td>Predictor Covariate</td><td>Description</td></tr></thead><tbody><tr><td>week</td><td>Time point (weekly ID)</td></tr><tr><td>driving</td><td>Vehicle increase rate (provided by Apple Inc.; compared with 13 January 2020)</td></tr><tr><td>transit</td><td>Public transport increase rate (provided by Apple Inc.; compared with 13 January 2020)</td></tr><tr><td>walking</td><td>Pedestrian increase rate (provided by Apple Inc.; compared with 13 January 2020)</td></tr><tr><td>num_beds</td><td>Number of available beds in hospitals in Tokyo (provided by Ministry of Health, Labour and Welfare)</td></tr><tr><td>vaccine</td><td>Number of people vaccinated in Tokyo (provided by Ministry of Health, Labour and Welfare)</td></tr></tbody></table>

These new features provide valuable information for investigating the dynamical system's behavior across various population segments. This facilitates the identification of similarities among segments for which the system exhibits identical behavior. Various classification techniques have been explored to achieve this functionality. More details, illustrative figures, and discussion can be found below on the experiments performed on the actual COVID-19 data.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** For each population segment, the link between control terms $\hat{B}_i, \hat{u}^*$, and $C^{(1)}, \dots, C^{(J)}$, which describe the characteristic variables (covariates, predictors) of the population segment, was fitted to a GAM. The z-score output of the GAM was used for classification to identify how these predictors contributed to the control of the system.</figcaption></figure>

### 3. Experiments and Application to Real Data

#### 3.1. *Data Description*

We used two sets of observed COVID-19 data: (i) data from five prefectures in Japan provided by JX PRESS Corporation and (ii) data from nine countries sourced from a "world-wide epidemiological database for COVID-19" [12] and the "COVID-19 data hub" [13]. Using our hidden controller framework, we estimate the unknown parameters $A, B$, and the gain matrix $K$ of the closed-loop system (9) to capture the dynamic behavior of COVID-19.

The prefectural data used in this study comprised the proportion of weekly infected individuals in prefecture $i$ ($\pi_i^1(k)$), the proportion of weekly recovered individuals ($\pi_i^2(k)$), the proportion of weekly deceased individuals ($\pi_i^3(k)$), and the proportion of all other
people who never had contact with the virus ($\pi_i^4(k)$(= 1 - $\sum_{n=1}^3 \pi_i^n(k)$)). The proportions were estimated from weekly counts of numbers reported daily for 32 weeks from 1 April 2021 to 10 November 2021, coinciding with the vaccination campaign in Japan. The five prefectures ($i = 1, \dots, 5$) comprise Tokyo, Osaka, Hokkaido, Fukuoka, and Okinawa, respectively. These prefectures are representative of Japanese populations for which there were a sufficient number of cases for statistical analysis. For each prefecture, the data are represented as a 32 × 3 matrix, with the weekly observations in the columns. The candidate weekly predictor covariates listed in Table 1 are for the period two weeks prior. Factors such as increases in driving and public transport (“transit”) in Tokyo, which are considered to affect the rate of new infections, are included. We used a two-week delay between measurement and onset in accordance with previous findings [14].

We selectively used data from nine countries, having complete data in our experiments: Australia (AUS), Brazil (BRA), Chile (CHL), Colombia (COL), Czech Republic (CZE), Germany (DEU), Lithuania (LTU), South Africa (ZAF), and Japan (JPN). The proportion of weekly infected individuals in country $i$ ($\pi_i^1(k)$), the proportion of weekly recovered individuals ($\pi_i^2(k)$), and the proportion of weekly deceased individuals ($\pi_i^3(k)$) were calculated from the data on “confirmed” (cumulative number of confirmed cases), “deaths” (cumulative number of deaths), and “recovered” (cumulative number of patients released from hospitals or reported recovered) [12,13]. For the sake of comparison, the period was the same as for dataset (i): the 32 weeks from 1 April 2021 to 10 November 2021. The candidate predictor covariates for these nine countries are shown in Table 2. Four of them represent governmental policy measures: government_response_index, stringency_index, containment_health_index, and economic_support_index. Since these policy measures are correlated and the ranges of the index scores vary between countries, we first categorized the policy measures into three levels—low (L: 0–50 percentile), middle (M: 50–90 percentile), and high (H: 90–100 percentile), using the Japanese data as a reference. We then applied multiple correspondence analysis and hierarchical clustering on principal components to the categorized data and obtained three groups of policy measures, as shown in Table 3. We again used a two-week delay between measurement and onset in accordance with the previous findings [14].

<table><caption>**Table 2.** Candidate predictor covariates for nine countries [12,13].</caption><thead><tr><td>Predictor Covariate</td><td>Description</td></tr></thead><tbody><tr><td>week</td><td>Time point (weekly ID)</td></tr><tr><td>driving</td><td>Vehicle increase rate (provided by Apple Inc.; compared with 13 January 2020)</td></tr><tr><td>transit</td><td>Public transport increase rate (provided by Apple Inc.; compared with 13 January 2020)</td></tr><tr><td>walking</td><td>Pedestrian increase rate (provided by Apple Inc.; compared with 13 January 2020)</td></tr><tr><td>hosp</td><td>Number of hospitalized patients</td></tr><tr><td>people_fully_vaccinated</td><td>Number of people who received all doses prescribed by vaccination protocol</td></tr><tr><td>government_response_index</td><td>Index of how government response varied over all indicators in Oxford COVID-19 Government Response Tracker database [15], becoming stronger or weaker over the course of the outbreak</td></tr><tr><td>stringency_index</td><td>Index of strictness of ‘lockdown’-related policies that restricted people’s behavior</td></tr><tr><td>containment_health_index</td><td>Index combining ‘lockdown’ restrictions and closures with measures such as testing policy, contact tracing, short term investment in healthcare, and investment in vaccines</td></tr><tr><td>economic_support_index</td><td>Index representing such factors as income support and debt relief</td></tr></tbody></table>
<table><caption>**Table 3.** Main categorical levels consisting of three groups of policy measures.</caption><thead><tr><th>Index</th><th>Government<br>_Response</th><th>Stringency</th><th>Containment<br>_Health</th><th>Economic<br>_Support</th></tr></thead><tbody><tr><td>policy 1</td><td>H</td><td>-</td><td>H</td><td>L</td></tr><tr><td>policy 2</td><td>M</td><td>H</td><td>-</td><td>M</td></tr><tr><td>policy 3</td><td>L</td><td>-</td><td>L</td><td>-</td></tr></tbody></table>

## 3.2. *Experimental Conditions*

To calculate parameters A, B from Equation (13) on the basis of the observed COVID-19 data for the five Japanese prefectures (i) and nine countries (ii), we first defined $Y = (a_{11}, a_{12}, \dots, a_{33}, b_{11}, b_{12}, \dots, b_{33})^T \in \mathbb{R}^{18=3*3+3*3}$, representing the unknown coefficients of parameters A, B. We then set the Hellinger distance (described in Section 2.3) to $d(x_k^*, x_k)$, where the term $x_k^*$ represents the estimated COVID-19 dynamics, and the other term $x_k = (\pi_j^1(k), \pi_j^2(k), \pi_j^3(k))' \in \mathbb{R}^3$ represents the observed COVID-19 dynamics on the $k^{th}$ day.

The proportions $\Pi(k)$ were calculated using the total population of Japan and the population ratios of each prefecture as of January 2021. We used the `fmincon` programming solver in MATLAB and selected the interior-point algorithm, for which the parameters of the termination tolerances on first-order optimality, function value, and input matrices were set to $10^{-4}$. The function optimized in `fmincon` was estimated for each window with a length of 3 using the $\ell_2$-norm. We used `mgcv` in R [10,11] to interpret the control parameters (Section 2.6).

The code and data pertaining to the nine countries are accessible via our GitHub repository: https://github.com/ismstat/hidden-controller-framework.git, accessed on 22 August 2023.

# 4. Results and Discussion

## 4.1. *Evaluating Predictive Power of Hidden Controller Framework*

Figure 5 shows comparisons between our hidden controller framework's prediction ($x_k^*$) and the observed $x_k$ for Australia, Colombia, Germany, and Japan based on the data for nine countries (data set (ii)). The predictions are reasonably accurate, suggesting that model parameters A, B were well-estimated. Similar results were observed for the other five countries, as detailed in Appendix A. These results underscore the predictive capacity of our hidden controller framework. Furthermore, control feedback $u_k$ can be effectively used in subsequent data analyses.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison between predictions made by our hidden controller framework (green dashed lines) and observed data (red solid lines) for (a) Australia, (b) Colombia, (c) Germany, and (d) Japan.</figcaption></figure>
## 4.2. Predictor Covariates Aimed at Halving the Number of Infected and Deceased Individuals

To analyze which predictor covariates effectively reduce the number of infected and deceased individuals, we applied our hidden controller framework (Section 2.3) and a multivariate GAM (Section 2.6) to a dataset in which the number of infected and deceased individuals was halved (dataset (i)). Figure 6 shows heatmaps depicting the z-scores of the predictor covariates associated with the numbers of (a) infected and (b) deceased individuals. Since z-scores approximately follow a Student's t-distribution and transition to a normal distribution when the degree of freedom becomes large, we calculated the z-scores at a significance level of 0.05, assuming a normal distribution. We color-coded the heatmaps blue for smaller z-scores, representing positive effects in decreasing the number of infected and deceased individuals, and red for larger z-scores, representing negative effects that increase these numbers.

We roughly categorized the predictor covariates into two groups, F1 and F2. F1 includes “week”, “transit”, “walking”, and “num_beds”, which generally had negative effects. In contrast, the predictor covariates in F2 (“driving” and “vaccination”) had positive effects. These results suggest that, in Japan, traveling by automobile and getting vaccinated effectively reduced the number of infected and deceased individuals.

We can classify the prefectures into two groups: the major metropolises of Tokyo, Osaka, and Fukuoka and the tourist cities in Okinawa and Hokkaido. The data suggest that vaccination in metropolitan areas is highly effective.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Heatmaps depicting z-scores of predictor covariates associated with numbers of (a) infected and (b) deceased individuals, with a goal of reducing the counts in both categories by half.</figcaption></figure>

To further investigate the positive and negative effects on reducing the number of infected and deceased individuals, we analyzed the difference in z-scores for predictor covariates associated with both categories, with the goal of halving each compared with the observed numbers. Figure 7 presents heatmaps illustrating these differences in z-scores. For the infected category, vaccination was highly effective in the metropolises of Tokyo, Osaka, and Fukuoka. For the deceased category, vaccination was particularly effective in Tokyo.

## 4.3. Differences in Predictor Covariates Across Countries

Figure 8 presents a heatmap of z-scores for predictor covariates associated with the number of infected, deceased, and recovered individuals across the nine countries. As an illustration of the classification of characteristics based on control (Section 2.6), the countries can be broadly categorized into four clusters: F1 includes South Africa, Chile, Brazil, and Colombia; F2 is comprised solely of Japan; F3 encompasses Australia and Lithuania; and F4 contains the Czech Republic and Germany. For group F4, three policy measures in Table 3 (P1 as indicated in the figure) and vaccination (P2) seem to effectively control the
number of infected, deceased, and recovered individuals. However, for countries in group F1, the opposite effects were observed.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Heatmaps illustrating difference in z-scores for predictor covariates associated with numbers of (a) infected and (b) deceased individuals, aimed at reducing each category by half compared with the observed numbers.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Heatmap of z-scores for predictor covariates (cov) associated with the number of infected ([cov].1), deceased ([cov].2), and recovered individuals ([cov].3) for nine countries.</figcaption></figure>

### 4.4. Clustering of Countries on the Basis of Policy Measures

Figure 9 illustrates the clustering of countries into four groups on the basis of the three categories of policy measures outlined in Table 3. These groups align with the ones determined by the z-scores in the previous section. This underscores the critical role of policy measures in controlling the number of infected, deceased, and recovered individuals.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Dendrogram illustrating clustering of countries on the basis of policy measures.</figcaption></figure>

**5. Conclusions**

We have presented a methodology that enables the correlation of various population characteristics with the evolution of a dynamically controlled system. We assume that the observed dynamical system is the result of the optimal control of a Kalman-type model, with no direct access to the various stages of this control. The choice of a Kalman-type model was made for simplicity, but this concept can be readily generalized to other models. We assume knowledge of only the model's structure, not its parameters. Our approach, therefore, is to estimate what the controller has carried out to arrive at the observed data; the estimates of the control parameters are made a posteriori. Our innovation resides in providing a data-driven solution, obviating the need for complex modeling. Furthermore, a distinctive facet of our work lies in treating the evolutionary objects of the dynamical system as probability distributions. We delve into exploring how these distributions have evolved under various hidden controller settings. We have provided an algorithm for estimating the various parameters that control the evolution of these distributions, which are based solely on a posteriori observations.

Once the control parameters were obtained, we used a multivariate generalized additive model (GAM) to explore the potential relationships between the part that controls the system's evolution and the various characteristics of the studied population. We applied this methodology to a COVID-19 database for five prefectures in Japan, considering various explanatory variables. We newly discovered that for the infected category, vaccination was highly effective in the metropolises of Tokyo, Osaka, and Fukuoka; while for the deceased category, vaccination exhibited particular effectiveness in Tokyo. This approach was generalized to nine countries in order to explore the effects of the various policies they follow. Our findings illustrate that countries can be broadly grouped into four clusters, with each group demonstrating varied effectiveness of policy measures and vaccination efforts.

From the GAM, we extracted indicators to quantify the degrees of association between the different explanatory variables. These new indicators enable the analysis and visualization of groups with similar and different evolutionary profiles. The aim is to explain the reasons for these differences. This work provides a novel perspective on understanding and managing complex evolutionary processes.

**Author Contributions:** Conceptualization, N.A.; Methodology, T.M. and N.A.; Software, T.M.; Validation, T.M. and N.A.; Formal analysis, T.M.; Investigation, T.M. and N.A.; Data curation, D.M.; Writing—original draft, T.M. and N.A.; Writing—review & editing, D.M.; Visualization, D.M.

**Funding:** This research received no external funding.

**Acknowledgments:** This research was supported by a "Strategic Programs" grant from ROIS (Research Organization of Information and Systems), Japan.

**Conflicts of Interest:** The authors declare no conflict of interest.
## Appendix A. Additional Results for Predictive Power of Hidden Controller Framework

Figure A1 presents a comparison between the predictions made by our hidden controller framework and the observed data for five countries. Consistent with the results in Section 4.1, our hidden controller framework performed well.

<figure><img src="image_13.png" /><figcaption>**Figure A1.** Comparison between the predictions made by our hidden controller framework (green dashed lines) and the observed data (red solid lines) for the remaining five countries.</figcaption></figure>

## References

1. Brogan, W.L. *Modern Control Theory*; Pearson Education India: Chennai, India, 1991.

2. Bryson, A.E.; Ho, Y.C. *Applied Optimal Control: Optimization, Estimation, and Control*; Routledge: London, UK, 2018.

3. Chen, C.T. *Linear System Theory and Design*; Saunders College Publishing: Philadelphia, PA, USA, 1984.

4. Fadali, M.S.; Visioli, A. *Digital Control Engineering: Analysis and Design*; Academic Press: Cambridge, MA, USA, 2012.

5. Bittanti, S.; Laub, A.J.; Willems, J.C. *The Riccati Equation*; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2012.

6. de Gouvêa, M.; Odloak, D. Dealing with Inconsistent Quadratic Programs in a SQP Based Algorithm. *Braz. J. Chem. Eng.* **1997**, *14*, 63–80. [CrossRef]

7. Wright, S.; Nocedal, J. *Numerical Optimization*; Springer Science: New York, NY, USA, 1999; Volume 35, p. 7.

8. Fletcher, R. *Practical Methods of Optimization*; John Wiley & Sons: Hoboken, NJ, USA, 2013.

9. Schittkowski, K.; Zillober, C. Nonlinear programming: algorithms, software, and applications. In *Proceedings of the IFIP Conference on System Modeling and Optimization*, Sophia Antipolis, France, 21–25 July 2003; pp. 73–107.

10. Wood, S.N. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. *J. R. Stat. Soc. (B)* **2011**, *73*, 3–36. [CrossRef]

11. Wood, S.N.; Pya, N.; Säfken, B. Smoothing parameter and model selection for general smooth models (with discussion). *J. Am. Stat. Assoc.* **2016**, *111*, 1548–1575. [CrossRef]

12. Guidotti, E. A worldwide epidemiological database for COVID-19 at fine-grained spatial resolution. *Sci. Data* **2022**, *9*, 112. [CrossRef] [PubMed]

13. Guidotti, E.; Ardia, D. COVID-19 Data Hub. *J. Open Source Softw.* **2020**, *5*, 2376. [CrossRef]

14. Matsui, T.; Azzaoui, N.; Murakami, D. Analysis of COVID-19 evolution based on testing closeness of sequential data. *Jpn. J. Stat. Data Sci.* **2022**, *5*, 321–338. [CrossRef] [PubMed]

15. Hale, T.; Angrist, N.; Goldszmidt, R.; Kira, B.; Petherick, A.; Phillips, T.; Webster, S.; Cameron-Blake, E.; Hallas, L.; Majumdar, S.; et al. A global panel database of pandemic policies (Oxford COVID-19 Government Response Tracker). *Nat. Hum. Behav.* **2021**, *5*, 529–538. [CrossRef] [PubMed]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.