1 Introduction
The modeling of time dynamical systems is of interest in multiple scientific fields. While ordinary differential equations (ODE) have a long history, interest in delay differential equations (DDE) is more recent, albeit both have been extensively studied
(Coddington and Levinson, 1955; Weiss, 1967; Ablowitz and Ladik, 1975; Driver, 2012; Asl and Ulsoy, 2003). A DDE is a natural extension of an ODE when observed processes have an aftereffect. Unlike the situation for ODEs, the solution of a DDE depends not only on the initial state but on the entire history of the process in a time interval of length equal to the delay prior to the initial time point.Moving beyond the basic notion of a deterministic ODE, random differential equations (RDE) (Strand, 1970; Soong, 1973; Cortés et al., 2007; Neckel and Rupp, 2013) and stochastic differential equations (SDE) (Arnold, 1974) are used to accommodate probabilistic uncertainty in temporal stochastic processes. In this paper we focus on RDEs, where the random effects are manifested in the model parameters. These include coefficients, initial conditions and forcing terms that are typically smooth in time, leading to differentiable sample path solutions of the RDE, in contrast to the situation for SDEs, which have a nondifferentiable forcing term, typically a Wiener process, and therefore require stochastic calculus (Itô, 1951).
Extensive developments in the theory of the forward problem of obtaining solutions for RDEs stand in contrast with statistical approaches, which include dataoriented methodology for the inverse problem, that is, learning the nature of the differential equation from data. This approach has been referred to as empirical dynamics (Müller and Yao, 2010). The starting point is a sample of random trajectories that are viewed as independent and identically distributed (i.i.d.) realizations of an underlying smooth (continuously differentiable) stochastic process. For a smooth stochastic process , there always exists a function such that
(1)  
(2) 
with almost surely. This forms the basis of empirical dynamics and has led to both parametric and nonparametric modeling approaches for (Zhu et al., 2011; Verzelen et al., 2012), which can be characterized as dynamics learning from data. Due to the minimal assumptions, dynamics learning covers a large number of specific ODEs that do not need to be a priori specified, however it requires that a sample repeated realizations of the underlying stochastic process is available. Such samples form the backbone of functional data analysis (Ramsay and Silverman, 2005; Wang et al., 2016). A pertinent example is provided by COVID19 caseload trajectories, where one observes samples of such trajectories across geographic units such as states or countries (Carroll et al., 2020a).
Related but quite distinct statistical inference concerns the problem of identifying the parameters of an a priori specified dynamic system (Brunel, 2008; Li et al., 2005; Chen and Wu, 2008), usually just from one or very few realizations. Here the trajectories are considered nonstochastic but are typically assumed to have been measured with noise. This leads to a curve fitting problem that in simple cases can be addressed with nonlinear least squares, but for which also nonparametric and semiparametric statistical methods have been employed (Liang and Wu, 2008; Paul et al., 2011). A key difference is that in dynamics learning one requires data that can be considered as an independent and identically distributed sample of realizations of an underlying smooth stochastic process , while the parameter identification problem is often addressed given (noisy) data from one trajectory. The modeling approach that we introduce here follows the paradigm of dynamics learning, albeit in a somewhat more structured framework than empirical dynamics.
The general form of DDEs for vector functions is
where is the function value or state at time and can be either a vector of the states evaluated at discrete time delays , so that , or alternatively an integral of the trajectory over a past period, representing a continuum of delays often called distributed delays (Bellen and Zennaro, 2013; Jacovitti and Scarano, 1993; Bjorklund and Ljung, 2003; Elnaggar et al., 1989; Mehrkanoon et al., 2014; Caraballo et al., 2018; Beretta et al., 2001; Gurney et al., 1980; Aziz and Amin, 2016). Results on existence and uniqueness of random differential equations with delay (RDED) are relatively recent, where Calatayud et al. (2019) studied conditions for existence and uniqueness of the solutions of a RDED and Cortés and Jornet (2020) the specific case of a linear RDED with a forcing term.
In this paper we propose dynamics learning from a sample of multivariate functional data, where the model generating the observed derivative trajectories is assumed to be a RDE with a delay component. The randomness in the DE is included in the forcing term, part of which is explained by additional covariates. These are stochastic processes with individual delay components, which in the motivating COVID19 application correspond to trajectories of mobility and economic activity. The unexplained remainder is a drift process, which appears as in equation (2).
We utilize tools from the theory of RDE to provide a foundation for the proposed population models and establish that all the models described in Section 2 correspond to RDED with unique solutions, which depend on the model parameters. The goal is then to learn these unknown parameters, which govern the dynamics encapsulated in the population RDED, by pooling information across the sample of observed trajectories. The presence of a large number of covariates and delay components to choose from gives rise to the challenge of model selection. To address this we propose an initial pruning step as described in Section 4.2 followed by a backfitting step to optimize over a set of delay components as outlined in Section 4.3. The proposed approach differs from existing optimization strategies (Wang and Cao, 2012; Zhou, 2016; Mehrkanoon et al., 2014; Wang and Wang, 2020) and statistical methodology (Jarne et al., 2017)
for parameter estimation that has been previously deployed to identify a DDE based on noisy data observed for a single nonstochastic trajectory. To our knowledge, dynamics learning where one has samples of stochastic processes has so far not been explored for the case of an underlying RDED, even for the case of a onedimensional stochastic process. We furthermore demonstrate in this paper that such an approach is well suited for modeling the growth rates of COVID19.
We consider both discrete and distributed delay models and establish existence and uniqueness of the solutions in the sense (Calatayud et al., 2019) in Section 3. For the case of discrete delays we harness functional concurrent regression models (Cai et al., 2000; Huang et al., 2004; Şentürk and Müller, 2008b; Şentürk and Nguyen, 2011; Huang et al., 2004) and for distributed delays the functional historyindex model (Malfait and Ramsay, 2003; Şentürk and Müller, 2010), which incorporates a range of recent past values of the process. For dynamics learning of RDEDs, we adopt a two stage procedure, where we first utilize functional linear regression (Cardot et al., 1999; Yao et al., 2005; Morris, 2015) with history index to learn the distributed delay, where the regression parameter function then corresponds to a history index function for the process of interest. In a second step the resulting linear predictor, which is the inner product of the history index function and the predictor process of interest, is used as a predictor along with additional covariate processes with covariatespecific delays in a concurrent model to fit the derivative process.
We apply this model to predict the timedynamic growth rate of COVID19 for individual states in the United States, where the sample of all states and their case trajectories provides the sample of trajectories that is the starting point for the proposed dynamics learning. Predictor processes that we consider include the cumulative case (caseload) process, daily economic indicators and changes in mobility patterns. Modeling the time evolution dynamics of the COVID19 pandemic is of great importance to understand and interpret the underlying associations as well as for deploying resources and formulating policies in the face of great uncertainty (Brett and Rohani, 2020; Bertozzi et al., 2020; Hao and Wang, 2021). The proposed methodology is also of interest to assess the dynamics of many other empirically observed complex multivariate stochastic processes for which samples of observed trajectories are available. We show by means of leaveoneout predictions that by employing dynamics learning for the proposed RDED model with delay components one obtains considerably more accurate timedynamic growth rate predictions for COVID19 caseload curves compared to models without delays, demonstrating the importance of considering the inclusion of lags when modeling the growth rate of COVID19.
2 Proposed Models
Let denote a multivariate stochastic process where is a continuously differentiable process of interest, is a vector function of additional covariates, and is a time window of interest. Consider a RDED with discrete delays,
(3) 
where corresponds to the initial condition stochastic process, the are discrete delays, and are smooth functions whose regularity will be specified below in Section 3. In the above, is a random drift process that is independent of .
While the inclusion of discrete delays is an extension of the classical functional concurrent regression model (Ramsay and Silverman, 2005), the functional linear regression model with history index (Malfait and Ramsay, 2003; Şentürk and Müller, 2010) might be a better choice when the derivative trajectories depend not only on the predictor process value at a single past instant but on the entire continuum in the recent past. We model these distributed delays as
(4) 
where again is an initial condition process. For the purpose of illustration and technical derivations, we assume that is a univariate process in (2); the corresponding multivariate generalization is straightforward.
For COVID19 modeling a hybrid model that combines distributed delay in with discrete delays for all other predictors turned out to be particularly suitable. In this practically relevant variant of a RDED one postulates that for all ,
(5) 
with an initial condition
We remark here that the RDED models proposed above are linear in the model parameters, although general RDED models determined by nonlinear functions are possible and might provide greater modeling flexibility, at the cost of additional modeling complexity. For example, it is not difficult to develop population models for quadratic and polynomial versions, in analogy to functional polynomial models (Yao and Müller, 2010). Modeling with linear RDEDs has the advantages that such models are easy to apply and implement, come with excellent interpretability of the model parameters, and exhibit good empirical performance in our application. For the rest of the manuscript we therefore only consider linear RDEDs and leave the development application of nonlinear RDEDs for future research.
3 Existence and Uniqueness of Solutions
3.1 Functional concurrent model with discrete delays
Here we discuss the existence and uniqueness of the solutions of the RDEDs corresponding to the concurrent model (2) and the historyindex functional linear model in (2), as well as the hybrid model in (5) used in the data application.
Consider a complete probability space
, where is the algebra on andis a probability measure. Denote the space of random variables with finite
moment by
, that is, for if . A sequence of realvalued random variables is said to converge to a random variable in the moment if as . We define a process to be moment continuous if for any The notions of moment differentiability and moment Riemann integrability are defined similarly (see Soong (1973) Chapter , pages and ).We say that a stochastic process is a solution of an RDED in the sense on the interval if is moment differentiable on , moment continuous on , and satisfies the RDED including the corresponding initial condition (Calatayud et al., 2019). We require the following assumptions on the coefficient functions and predictor processes in model (2).

) and are continuous on .

is moment continuous on .

The predictor processes satisfy and are moment continuous on , for
Calatayud et al. (2019) showed the existence and uniqueness of the solution of the following general form of a random delay differential equation with discrete delay , given by
(6) 
and established the following result.
Theorem of Calatayud et al. (2019). If satisfies the Lipschitz condition: for and a with then the random delay differential equation (3.1) has a unique solution in the sense.
The concurrent model in (2) is a special case of the random delay differential equation (3.1) with a specific choice of . A key step is to show that the function defined in (2) is moment continuous. The existence and uniqueness of the solutions of (2) can then be obtained following similar arguments to those in the proof of Theorem 2.2 in Calatayud et al. (2019). We summarize this result in Proposition 1. All proofs are in Appendix A.1.
Proposition 0.
Remark 0.
Proposition 1 requires the moment continuity of on for the uniqueness of the solution in the sense. If almost surely, then there exists a version of the solution which solves the RDED in the sample path sense, see Theorem 2.3 in Calatayud et al. (2019) and Caraballo et al. (2019). Note that in (2), we do not include the process itself as one of the predictors; included is only the delayed version , where . Hence solutions in the sample path sense do not require to be differentiable almost surely. For , see Remark 4.
Remark 0.
The explicit sample path solution in the special case of a linear RDED has been derived in Cortés and Jornet (2020); Caraballo et al. (2019); Calatayud et al. (2019). The linear RDED often includes random coefficients and a random timevarying forcing term. In (2), the terms involving the predictor processes , the drift process , and the intercept coefficient function form the constituents of the random forcing term; the varying coefficient functions are nonrandom. Using the method of steps as outlined in the proof of Theorem 2.3 in Calatayud et al. (2019), we can provide a solution to (2) in the sample path sense as follows. For , (2) is equivalent to an ordinary RDE given by
(7) 
where is defined as in (2). Then, the solution of (3) is given by
Next, for ,
(8) 
Solving (3) leads to
Repeating this argument until the entire domain is covered one obtains the sample path solution. Due to the existence and uniqueness of the solution of (2) in the sense and by Remark 2, the sample path solution obtained as above is unique.
Remark 0.
In applications, the delay parameter can take the value zero and when in (2) one arrives at an ordinary RDE without delay. For this case, existence and uniqueness of solutions has been well studied (see, e.g., Theorem 5.1.1 in Soong (1973) and Strand (1970)). The sample path solution in this case can be obtained via the integration factor method and is given by
where with defined as in (2). Here consists of the intercept term , the predictor processes along with the corresponding coefficient functions , and the drift process .
3.2 Functional linear model with history index for distributed delay
We start with the simplest history index model with only one predictor process , characterized as a first order RDED described in (2). We prove the existence of a unique solution of (2) by adopting similar arguments as those behind Proposition 2.1 and Theorem 2.2 of Calatayud et al. (2019). The extension to multivariate predictor processes is straightforward. Observe that one can rewrite model (2) as
(9) 
where with the initial condition
Assume that is moment continuous on and define such that for a moment continuous process and ,
(10) 
Here is a random functional whose first argument is the trajectory and whose second argument is . The following proposition with proof in Appendix A.2 provides the continuity of in in the sense.
Proposition 0.
Suppose is moment continuous for . Then as defined in (10) is continuous in in the sense on the domain .
For the random process on and the time variable , the history index model (9) can be expressed as
(11) 
For the existence of a unique solution of (11) we need the following assumptions, in addition to a c:

is continuous on and is continuous on .

on and is moment continuous on .
The following result characterizes the solution of (11), where the integral in condition b of the following Proposition is defined in the moment Riemann integral sense (for definition of moment Riemann integrability see Soong, 1973, Chapter page ).
Proposition 0.
The process is a solution of (11) if and only if for

is moment continuous on

for each

for each .
The following result establishes existence and uniqueness of a solution for (11). The proof extends the classical Picard theorem for deterministic ODE to RDED in the sense, via the Banach fixedpoint theorem, following a similar line of arguments as in Calatayud et al. (2019) (see Appendix A.2).
Theorem 7.
For the existence and uniqueness of the solution of the hybrid model in (5), which is the most promising of the models considered for modeling the growth rate of COVID19, the following corollary applies.
4 Learning Dynamics from Samples of Multivariate Stochastic Processes and Practical Considerations for COVID19 Caseload Modeling
4.1 Obtaining derivatives
In datadriven differential equation modeling, given a sample of observed processes, a necessary first step is the estimation of derivatives. In applications such as the COVID19 case trajectories the available data are usually noisecontaminated, and therefore some care is needed to obtain viable derivative estimates. The situation for the observed data for one of the trajectories is reflected by the nonparametric regression model
(12) 
where are the grid points where the process is measured and
are the measurement errors, which typically are considered to be independent mean zero random variables with finite variance.
In our data application, the available data correspond to the cumulative case counts of COVID 19, which we consider to be noisy realizations of a smooth underlying function that is observed on a daily grid. For the practical estimation of the derivatives one has various choices that all require a tuning parameter. Since in preliminary studies local polynomial fitting yielded the smallest leaveoneout prediction error relative to the difference quotients, which serve as nearly unbiased targets, we opted to obtain derivatives by applying local polynomial regression (Müller et al., 1987; Fan and Gijbels, 1996) as per Section A.3 of the Appendix, where further details are provided.
4.2 Target model and initial lag and covariate process selection
For modeling the COVID19 caseload process for the states of the U.S. we adopt the RDED model (5),
where the time delay of the effect of the process at time is distributed on the interval , and are discrete time delays or time lags for the covariate or predictor processes . Our starting point is that the data represent independent realizations of the stochastic process , where each realization (corresponding to a state) is indexed by and is the number of states included in the analysis.
Given , to estimate the history index weight function in the historical model without covariates, we preform a scalartofunction linear regression for each grid point , i.e.,
(13) 
This can be implemented using function FLM in the R package fdapace (Carroll et al., 2020b). With estimates of for each in hand, the estimate of the history index weight function, for all , is obtained using a two dimensional kernel smoother.
For each of the covariates , we conduct an initial lag selection based on functional varyingcoefficient regression including the corresponding variable as a single predictor, considering the models
(14) 
where . Given a lag , for each subject and time point , we obtain the leaveoneout prediction
where and are estimated by fitting model (14
) by ordinary least squares, excluding the data from subject
. The lags are chosen to minimize the leaveoneout prediction error across subjects,(15) 
With initial choices of the lags and the estimate obtained from model (13), for each , variable selection is then performed using LASSO (Tibshirani, 1996),
where . This can be implemented using the R package ncvreg (Breheny and Huang, 2011; Breheny, 2020), where the tuning parameter is chosen by leaveoneout crossvalidation. Considering all time points simultaneously, we obtain proportions of the number of time points where each variable is selected, i.e.,
Subsequently, selection among the timevarying predictors other than the process is performed by applying a threshold on these proportions . This threshold is chosen by leaveoneout crossvalidation,
(16) 
4.3 Lag selection by backfitting
Let denote the set of selected predictor processes obtained using the method described in Section 4.2. Updated lag parameters are then chosen for the variables in through an iterative backfitting algorithm such that the total leaveoneout prediction error is minimized. Given the iteration of the lag vector, one obtains the update for the lag by minimizing the leaveoneout integrated mean squared prediction error, i.e.,
where for
Here the functions represent the fitted coefficients corresponding to the iteration of the lag, with the observation held out for prediction. This process is iterated over until the lag vector converges. Lags are initialized at the vector , where are as described in Section 4.2 and the stopping rule is such that the algorithm terminates if no lags have changed after a complete cycle of updates for all predictors. After the final lags are estimated, we lightly smooth the estimated coefficients with local linear smoothing, which is a standard procedure in concurrent functional regression modeling (Fan and Zhang, 2000; Şentürk and Müller, 2008a). Specifically, the initial estimates are smoothed to obtain the effect functions presented in Figures 3 and 4, defined as where
(17) 
with representing the Epanechnikov kernel with a bandwidth of days.
5 Random Differential Equation with Delay Fitting and Prediction for COVID 19 Case Trajectories in the United States
5.1 Data Description
We obtained daily confirmed cases across states in the United States from the COVID19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University, which is publicly available at https://github.com/CSSEGISandData/COVID19 and was accessed on September , . The positive and negative COVID19 test results per day and state were obtained from the COVID Tracking Project and are publicly available at https://COVIDtracking.com/ (accessed on September , ). The former variable refers to the number of people with confirmed or probable case of COVID19 (see https://COVIDtracking.com/aboutdata/datadefinitions for more details) and the latter refers to the total number of people with a completed negative PCR test. We used effective positivity rates (EPR), defined as the ratio of effective positive results (positive test or probable case) to the effective total tests (positive or probable case plus negative test results).
We also obtained features from the Opportunity Insights Economic Tracker Data (Chetty et al., 2020), which is publicly available at https://github.com/OpportunityInsights/EconomicTracker (accessed on September , ). This database contains daily and state level information for several economic activity indicators such as credit and debit card expenditure across different types of activities, some of them stratified by zip codes with low/medium/high income level. We refer to Chetty et al. (2020) for more details. Additionally, we downloaded Google mobility data, publicly available at https://www.google.com/COVID19/mobility/ (accessed on September , ). These data include indicators of mobility pattern changes in different areas such as parks, residential locations, retail locations, among others; and information such as the percent change in revenue as well as total number of small businesses for high/low income consumers and across different economic activities; see Table 1 for a complete listing.
Item  Feature description 
1  cases per million 
2  ratio of effective positive results to effective total tests 
3  ratio of effective positive results to state population in 
4  individuals who are currently hospitalized with COVID19 
5  spending in all categories 
6  spending in arts, entertainment, and recreation categories 
7  spending in accommodation and food service categories 
8  spending in general merchandise stores and apparel and accessories categories 
9  spending in grocery and food store categories 
10  spending in health care and social assistance categories 
11  spending in transportation and warehousing categories 
12  spending among high (top quartile) income ZIP codes in all categories 
13  spending among low (bottom quartiles) income ZIP codes in all categories 
14  spending among middle (middle two quartiles) income ZIP codes in all categories 
15  time spent at retail and recreation locations 
16  time spent at grocery and pharmacy locations 
17  time spent at parks 
18  time spent at inside transit stations 
19  time spent at work places 
20  time spent at residential locations 
21  time spent outside of residential locations 
22  number of small businesses open 
23  number of small businesses open among high(top quartile) income ZIP codes 
24  number of small businesses open among low(bottom quartile) income ZIP codes 
25  number of small businesses open among middle(middle two quartiles) income ZIP codes 
26  number of small businesses open in transportation 
27  number of small businesses open in education and health services 
28  number of small businesses open in leisure and hospitality 
29  net revenue for small businesses 
30  net revenue for small businesses among high(top quartile) income ZIP codes 
31  net revenue for small businesses among low(bottom quartile) income ZIP codes 
32  net revenue for small businesses among middle(middle two quartiles) income ZIP codes 
33  net revenue for small businesses in transportation 
34  net revenue for small businesses in education and health services 
35  net revenue for small businesses in leisure and hospitality 
We briefly describe a few features below and refer to the Economic Tracker source (Chetty et al., 2020) for further details:

‘spend all’: Seasonally adjusted credit and debit card spending relative to the period Jan to Jan , , in all merchant categories, reported as a seven day moving average formed by using the values for the current and the previous six days.

‘spend tws’: Seasonally adjusted credit and debit card spending relative to the period Jan to Jan , , in the category transportation and warehousing, and also reported as a seven day moving average.

‘spend all inchigh’: Seasonally adjusted credit and debit card spending by consumers living in high median income (above quantile) zip codes, and relative to the period Jan to Jan , , in all merchant categories, again reported as a seven day moving average.

‘revenue all’: Seasonally adjusted percent change in net revenue for small businesses and indexed to the period Jan to Jan , , reported as a seven day moving average.

‘revenue ss70’: Seasonally adjusted percent change in net revenue for small businesses in leisure and hospitality and indexed to the period Jan to Jan , , reported as a seven day moving average.

‘merchant inchigh’: Seasonally adjusted percent change in the number of small businesses open in high median income (above quantile) zip codes and indexed to the period Jan to Jan , , reported as a seven day moving average.

‘gps parks‘: Corresponds to the time spent at parks relative to a baseline which is defined as the median value, for the corresponding day of the week, during the period Jan –Feb , .
We employ local linear regression to slightly smooth the feature trajectories; this also allows to impute missing observations for some states and days. For this step, a Gaussian kernel with bandwidth
days was employed and we utilized the R package fdapace version 0.5.5 for the computational implementation (Carroll et al., 2020b). The effective test positivity rate was smoothed analogously. We use local quadratic regression (Fan and Gijbels, 1996) to obtain the derivatives , implemented by using the R package locfit version 1.59.4 (Loader, 2020). For data analysis we consider states with population at least 1 million as the states with smaller populations had less COVID19 cases, and their inclusion would therefore increase the overall noise level in the data and negatively impact the analysis. This is a pragmatic choice and is not a limitation of the proposed method. In the end, we included the data from 44 states (excluding Alaska, Delaware, North and South Dakota, Vermont and Wyoming).As a result of the LASSO variable selection described in Section 4.3, which was implemented with the R package ncvreg version 3.12.0 (Breheny, 2020) with threshold in (16), the features utilized in the final model are as follows, where the selected lags are listed in Table 2.

, the total cases per million at time for state , .

, the relative time spent at park areas at time in state .

, the number of hospitalized people at time in state .

, the “positive testing rate” defined as the ratio of effective positive results (positive test or probable case) to the effective total tests (positive or probable case plus negative test results) at time for state .

, the credit/debit card relative spending in the arts, entertainment, and recreation category (see ‘spend aer’ predictor in Table 1) at time in state .
predictor  variable  lag  
park mobility activity  0  
number of people currently hospitalized  14  
effective test positivity rate  7  

3 
5.2 Model Fitting and Results
In order to study the delay time dynamics of COVID19 in the United States, we focus on the statespecific timevarying cumulative confirmed cases per million people, , for , which is referred to as the case process in the following. Here, the time domain that we consider starts on April 5 and ends on August 12 in 2020, where data are recorded once per day, on an equidistant grid for 130 days. We estimate the history coefficient function in the RDED model (5) as described in Section 4.2 with , and adopt the functional varyingcoefficient regression in model (5) with other variables as listed in Table 1. The estimated coefficient surface and some of its crosssections are displayed in Figure 1.
For each fixed , , the estimated history index in model (5) is found to change from positive near current time to negative with increasing lags. This means that the recent past of the process has a positive association with the current derivative while days further into the past are negatively associated. The right panel shows that the shape of changes from more or less quadratic to linear from the early days of COVID19 to recent times, indicating that in April 2020 the contrast between the effects of nearcurrent and past caseload on the current derivative was more pronounced than it is more recently.
5.3 Predictor process effects on growth rate
The results are visualized in Figure 2. States which are suffering from faster spread of the virus are characterized by higher growth rates. The historically weighted integrated case count per million has a positive effect on the growth rate. This reflects dynamic explosive behavior (Müller and Yao, 2010), so that states will tend to have caseloads that move even further away from the mean caseload per million across states as time moves forward, which could be either moving further above the average across states or further below it. Thus the positive effect of the history index reflects increasing variability in the caseloads across states. Potential reasons for this are the distinct policies that states implement in terms of mobility, business and social restrictions, as well as cultural acceptance of such restrictions. Spatial infectivity dynamics with waves of infections moving spatially across the U.S. also might play a major role.
In terms of the effects of covariate processes, our results indicate that higher park mobility during late spring and summer (May onwards) is associated with a decrease in growth rates. For each percent increase in park mobility, the growth rate decreased by as much as 10 cases per million per day, with the effect reaching its maximum in midsummer. This indicates negligible risk and potential benefits from outdoors recreation activities.
The other covariate processes that were selected by LASSO as predictors in the RDED model (5) are patient counts at hospitals, the testing positivity rate, and credit card spending on recreational activities, including dining at restaurants. All of these were found to have positive effects on the growth rate throughout nearly the entire time domain. These effects were included in model (5) with covariate processspecific lags, as per Table 2. The number currently hospitalized was found to be most predictive for the case growth rate when including a twoweek lag, whereas for test positivity rate a oneweek lag was found to be most predictive. This may suggest that individuals who are infected by individuals with less severe symptoms (i.e. a nonhospitalized COVIDpositive person) tend to learn of their disease one week after their infector did.
The twoweek difference for hospitalized individuals is likely a reflection of the time needed to develop serious symptoms which warrant admission for medical care. The effect of recreational spending on viral spread was much more immediate; the optimal lag in this case is only 3 days. This suggests that states with higher recreational economies during the pandemic see an almost immediate uptick in growth rate, with cases increasing by as much as 25 infections per million per day for each percent of additional recreational spending.
5.4 Predictive performance of the random differential equation with delay model
Figures 3 and 4 display the fitted growth rates provided by the RDED model (5), with covariate specific lags, and also the results of fitting a differential equation model where the covariates do not have individual lags, i.e. a concurrent RDE model given by (5) in which . Fitted curves are the result of leaveoneout predictions for each state, where the data from the leftout state were not used for fitting the model and then the state’s data are plugged into the fitted model to obtain the predicted case growth rate. These predictions are compared against the actual observed growth rates as obtained by local polynomial regression (see Section 4.1 and Section A.3 of the Appendix).
Visual inspection (as well as measures of goodnessoffit, such as integrated mean squared prediction error) indicates that the incorporation of delays with predictor specific lags in the model substantially improves the fits. This provides compelling evidence for modeling COVID19 growth rates with RDED that include covariatespecific lags. A vanilla RDE without predictor delays does not incorporate historical information beyond the previous day and struggles to explain deviations in phase (as in, e.g. Illinois). The problem of modeling phase variation is instead foisted upon the coefficient functions, which are illequipped to reflect it, and thus at times the straight RDE model without predictor delays falters in terms of amplitude modeling as well (e.g. for Colorado). In contrast, the proposed RDED incorporates past information through the included delays and thus reflects phase variability in a natural way, especially through the inclusion of distributed delays via the history index. Including a twoweek history of case counts through the history index function in model (5) substantially enhances the flexibility of the model. The estimates of coefficient functions are no longer uncontrollably altered by the effects of unaccounted lags, and as a result the fitted model performs much better.
5.5 Modelbenchmarked evaluation of performance
A comparison of observed vs. RDED modelpredicted growth rates allows one to assess the performance of states in terms of controlling spread. A state whose observed growth rate falls below the predicted curve corresponds to a negative residual and suggests betterthanexpected viral control, and vice versa for observed rates which fall above the predicted curve. Figures 5 and 6 display the residual curves for the fitted RDED model (5), where positive residuals correspond to the number of excess new cases per day relative to the predictions obtained from fitting model (5), as described above, including the twoweek history index and other covariate processes with their respective lags. These residuals reflect the unexplained residual stochastic process that is an integral and principally unpredictable part of the RDED model (5) on which our predictions are based.
Major outbreaks, which can be interpreted as an unforeseen amount of excess cases, are revealed through positive spikes in the residual plots. For example, New York’s residuals leap up to as high as 170 excess cases per million per day during a spike in Spring. Arizona’s summer outbreak is captured in the residuals as well, subsequently coursecorrecting with a streak of fewerthanexpected new cases around day 100. These spikes also suggest that in the two weeks prior to the uptick, there were substantially more cases in the state population than the numbers reported in the JHU dataset.
The consistency of a state’s handling of the virus can be seen through the volatility of the residual curve. States like Colorado and Oregon exhibit only minor deviations from zero, which suggests their new cases are quite predictable given a twoweek history of the state case counts and covariate processes. Some states, particularly in the Northeast, including New York, New Jersey, and Rhode Island, only exhibit high volatility during the early days of the pandemic before achieving a period of stability. On the other hand, states such as Alabama and Louisiana exhibit lesspredictable case counts throughout.
5.6 Computation times and code availability
The codes for reproducing the above analysis are available at https://osf.io/w48zd/?view_only=7cd3e2c686a74fd086e4eae8534214c0. The run times of the different steps in the analysis using R version 3.6.3 (20200229) running under CentOS Linux 7 (Core) on x86_64pclinuxgnu (64bit) platform are summarized in Table 3.
Step  Description  Number of Cores  Time in Seconds 
1  preprocessing  single core  9.475 
2  initial lag selection  8 cores in parallel  422.505 
3  history index function estimation  single core  3.414 
4  lasso  8 cores in parallel  31.316 
5  backfitting  single core  2565.715 
6  smoothing  single core  0.012 
Total  3032.437 
6 Discussion and concluding remarks
In this paper we propose dynamics learning exemplified by learning a RDED from a sample of observed trajectories, which is motivated by modeling COVID19 daily new cases from the previous behavior of the case trajectory itself through a distributed delay or history index component and additional covariate process components. RDEDs form the driving mechanism of real life processes that include a feedback loop into the future. Depending on the application and the exact nature of the feedback, one might incorporate discrete delays or distributed delays. While discrete delays borrow information from some isolated times in the past, a distributed delay takes into account the entire continuum in the recent past and turns out to be closely related to a functional history index model that has been considered in the statistics literature.
In the current framework, we have not included inference for the estimated model coefficient functions and the lags. It is left for future research to investigate the convergence of these estimated parameters to their population targets and to obtain guarantees on the convergence rates. Some words of caution are in order. Any method that includes a time delay component suffers from the limitation that a larger set of possible lags means one has to sacrifice a part of the data that are available for the analysis, due to the initialization. We have assumed that the trajectories corresponding to the different states are i.i.d, realizations, which may not hold if there is significant spatial correlation based on geographic and demographic similarities across the states. It would be interesting to study adaptations of the proposed methodology when the samples are correlated instead of being independent. From an application point of view, a different and perhaps more refined perspective could potentially be gained by applying the proposed method at a finer level to county level data.
Establishing existence and uniqueness of the solutions of RDEDs has found increasing interest in recent years (Calatayud et al., 2019; Cortés and Jornet, 2020). We contribute here to the forward problem of obtaining solutions of an RDED by targeting the inverse problem analytically and combine this with functional data analysis approaches to estimate the parameter functions from an observed sample of multivariate stochastic processes for both discrete and distributed delay setups. Our approach builds on and contributes to empirical dynamics and the literature on inferring dynamics and differential equations from functional data, and extends it to models that include time delays.
A key feature of the proposed approach is the statistical estimation of the model parameters in the population RDED, which has inbuilt model selection steps that aid in estimating the optimal lags and also to select the important predictors through the lasso pruning steps, with majority voting across all time points. This reduces the number of predictors in the model and thus its complexity before the optimization step of the lags that are associated with the predictors. To handle this complex lag optimization step, instead of a brute force search method, we adopted a backfitting algorithm, which significantly simplifies the optimization step and also leads to faster convergence. We provide a complete toolkit that starts with a sample of multivariate trajectories and extracts from the observed data an estimated RDED that best explains the dynamics inherent in the data, along with automatic variable selection and delay estimation.
The application of our method to the modeling of COVID 19 growth rates in the United States demonstrates excellent performance in terms of predictions of the trajectories. The predictions from the random differential equation model without predictor lags grossly overestimates the COVID19 growth rates for significant periods of time in New York, New Jersey, Pennsylvania, Massachusetts, Connecticut, Virginia, West Virginia, New Hampshire, California, Texas, Colorado, Oregon, Washington, Maine and Hawaii but underestimates the growth rates for Alabama, Arizona, Arkansas, Florida, Louisiana, Nebraska, Nevada, South Carolina, Tenessee and Wisconsin; see Figures 3 and 4. This highlights the importance of incorporating time lags in differential equation modeling of functional data, when one has a sample of realized stochastic processes, specially when the underlying processes tend to have an aftereffect.
Appendix
Appendix A.1 Proof of the existence of a solution for the model in (2)
Proof of continuity of .
The unique solution of the first order DDE in (2) can be shown to exist as a direct extension of Theorem 2.2 of Calatayud et al. (2019).
We first show that is moment continuous, using assumptions a  c in Section 3.1. Observe that
(18) 
The first and third terms in (A.1) converge to using assumptions a and b respectively. For the second and fourth terms of (A.1),
(19) 
By assumption a, is continuous on the compact domain , and hence uniformly bounded. Also, by assumption , in the moment. Since , the space of random variables with finite moments, is a Banach space, it follows that Thus (A.1) converges to . Next
Comments
There are no comments yet.