# Learning from limited temporal data: Dynamically sparse historical
functional linear models with applications to Earth science

Joseph<EMAIL_ADDRESS>[ Shizhe
<EMAIL_ADDRESS>[ Asad
<EMAIL_ADDRESS>[ Stefan
<EMAIL_ADDRESS>[ Jiguo
<EMAIL_ADDRESS>[ William J<EMAIL_ADDRESS>[
Nadja<EMAIL_ADDRESS>[ Ali A<EMAIL_ADDRESS>[ Department of Earth, Ocean and Atmospheric Sciences, University of British
Columbia presep=, ]e1,e3,e8  Department of Mathematics, University of British
Columbia presep=, ]e2  Department of Data Science, Norwegian University of
Life Sciences presep=, ]e4  Department of Statistics and Actuarial Science,
Simon Fraser University presep=, ]e5  Department of Statistics, University of
British Columbia presep=, ]e6  Norman B Keevil Institute of Mining
Engineering, University of British Columbia presep=, ]e7  Corresponding
author: Jiguo Cao<EMAIL_ADDRESS>
###### Abstract

Scientists and statisticians often want to learn about the complex
relationships that connect two variables that vary over time. Recent work on
sparse functional historical linear models confirms that they are promising
for this purpose, but several notable limitations exist. Most importantly,
previous works have imposed sparsity on the coefficient function, but have not
allowed the sparsity, hence lag, to vary with time. We simplify the framework
of sparse functional historical linear models by using a rectangular
coefficient structure along with Whittaker smoothing, then relax the previous
frameworks by estimating the dynamic time lag from a hierarchical coefficient
structure. We motivate our study by aiming to extract the physical rainfall-
runoff processes hidden within hydrological data. We show the promise and
accuracy of our method using four simulation studies, justified by two real
sets of hydrological data.

Functional data analysis ,

Hydrology,

Precipitation,

Rainfall-runoff relationships,

Streamflow,

Time lag,

Whittaker smoothing,

###### keywords:

, , , , , , and ,

##  1 Introduction

In many complex systems, the impact of a cause on some response of interest is
distributed across time. Exactly how this impact spreads and changes across
time is of great interest to scientists and policy-makers in a wide variety of
disciplines. This question not only relates to the magnitude of the effects,
but also to exploring how far back in time the cause has a significant impact
on the response. If a method for learning this impact function in a data-
driven manner existed, one would be able to extract useful knowledge from a
variety of time series data sources.

More specifically, this research is motivated by the desire to learn the data
generating processes behind hydrological temporal data (Botter, Bertuzzo and
Rinaldo, 2010). Learning the different data generating processes across a
variety of catchments around the world will allow hydrologists to infer the
geological or climatological factors that cause process variation (Schuite et
al., 2019). A greater understanding of process variation will enable better
predictions in regions where streamflow is not measured. Observational units
in hydrology called catchments, filter rainfall and snowmelt into streamflow
(Figure 1). This process may seem simple; however, slightly different
geologic, climatic, or topographic features may produce vastly different
filtering functions (Blöschl et al., 2013). The percent of precipitation that
is partitioned into streamflow versus water storage versus evaporation can
form wildly different ratios, and the time it takes for input water to become
streamflow can vary from minutes to years (Jasechko et al., 2016). Figure 1
illustrates these processes. After precipitation falls, it can either
transport out of the catchment via evapotranspiration, or it can be
transferred into flow. The depth to which water flows or gets stored in/on the
ground dictates the length of time between a precipitation event and a rise in
streamflow. Overland flow (OF) is transferred into streamflow in a few hours
or days. Shallow subsurface flow (SF) is transferred to streamflow in weeks or
months, while deep groundwater flow (GF) gets transferred in months or years.

Figure 1: A flow path diagram depicting how precipitation (red arrows) gets
partitioned into evapotranspiration (ET, yellow), or different types of flow
(green, blue, brown) on a hillslope. Overlandflow (OF, green), shallow
subsurface flow (SF, blue) and deep groundwater flow (GF, brown) eventually
reach the stream and become streamflow (light blue). The size of each arrow
depicts the relative influence of each component in a hypothetical catchment.

Figure 2 displays the daily streamflow and rainfall from the Koksilah River
(British Columbia) and the Withlacoochee River (Florida) for the year 1998,
illuminating the fact that the complex relationship between rainfall and
streamflow not only varies across different catchments, but also across
different seasons (Chiu and Bittler, 1969). At the Koksilah River, winter
rainfall does not get evaporated (days 1-50 and days 320-365), rather it
quickly transfers to streamflow. On the other hand, in summer and early fall
(days 100-310), any rainfall that may occur does not seem to result in notable
streamflow during these months (top of Figure 2). At the Withlacoochee River,
streamflow falls less rapidly during periods of little-to-no rainfall compared
to the Koksilah River (e.g., days 1-50 and days 80-120). Evaporation is likely
to be large in this catchment, therefore, streamflow may only increase in the
event of a large rainfall event as seen on day 50. In the summer and fall
(days 150-320), moderate to large rainfall events may not produce streamflow.

Figure 2: The daily streamflow (red) and rainfall (black) for the Koksilah
River (top) and the Withlacoochee River (bottom) for the year 1998.

Understanding these types of hydrological processes is vitally important for a
variety of reasons. For example, the mining sector, which is a significant
economic driver in many countries globally, is also highly vulnerable to local
hydrology. In Australia, flooding events in the Bowen Basin coal mining region
were attributed to billion-dollar production losses for the industry and led
to concerns about the impacts of saline mine discharge on surrounding
ecosystems (Vink and Robbins, 2012). Mining also alters its surrounding
landscape which can impact the water cycle. There are numerous examples of
mining regions that have perpetual water quality legacy issues, including
selenium pollution and acid rock drainage, which are mobilized due to local
hydrology (Kunz, 2020). Improved understanding of hydrological relationships
can serve an important role in mitigating such impacts. Likewise,
understanding the underlying hydrological structure of a landscape can serve
to identify regions that might be more hydrologically vulnerable to major land
use changes. For example, knowing that a catchment filters most rainfall
slowly through subsurface and deep flow pathways would indicate that the
stream is vulnerable to soil and deep rock contamination (Figure 1).

When predicting streamflow, it is common to assume an arbitrary fixed window
of time for which precipitation can affect streamflow, however, this is a
severe limitation of previous works (Janssen, Radić and Ameli, 2021), so
finding good constant or dynamic time lags is of interest to the hydrological
community. Further, learning these parameters at different catchments around
the world would allow hydrologists to learn what catchment factors have the
biggest influence on the shape of the impact function. Since climate and land
use changes will continue to occur, it would be of particular interest to
understand where these changes will have outsized impacts on local hydrology.
Before advancements in hydrology can be made in this direction, we must first
be able to accurately estimate the impact function of precipitation on
streamflow.

Xun, Guan and Cao (2022) considered the historical functional linear model,

$y_{i}(t)=\int_{0}^{t}\beta(s,t)x_{i}(s)\ ds+\epsilon_{i}(t),$ (1)

to study a constant time lag historical effect. Here $y_{i}(t)$ is the outcome
of interest at time $t$ of replicate series $i$, and $x_{i}(s)$ is a temporal
explanatory variable at the past time $s$ for the same replicate. The
bivariate coefficient function $\beta(s,t)$ represents the effect of the
explanatory variable at time $s$ on the outcome at the current time $t$, where
$0\leq s\leq t$. Note that for the entirety of this work, we ignore the
intercept term without loss of generality by removing the seasonal average
from $y_{i}(t)$ and $x_{i}(t)$. The residual term is denoted by
$\epsilon_{i}(t)$. To estimate $\beta(s,t)$, Xun, Guan and Cao (2022)
minimized an objective function containing three terms: (1) a least-squares
loss, (2) a nested group bridge penalty term, and (3) a smoothness penalty
term. These terms will be further discussed in Section 2, but in summary, the
two penalty terms allow the method to infer a sparse and smooth version of the
coefficient function, while the least-squares term pushes the optimization
towards a solution that fits the data well. More specifically, the group
bridge penalty forces some higher lag coefficients of $\beta(s,t)$ to zero,
thereby revealing the constant lag $\delta$, after which $\beta(s,t)$ is zero
and the covariate has no effect on the response. Sparse estimation of the
coefficient function is often of great interest, though it is rarely
considered (Kim, Şentürk and Li, 2011).

This very recent work by Xun, Guan and Cao (2022) builds upon many previous
methods including Malfait and Ramsay (2003) who originally developed the
functional historical linear model, and Harezlak et al. (2007) who proposed to
impose a discrete difference penalty on the coefficients similar to the
P-spline framework (Eilers, Marx and Durbán, 2015). In each of these works,
triangular finite elements with tent-like basis functions are considered to
model the coefficient function $\beta(s,t)$. The method introduced by Xun,
Guan and Cao (2022) is compelling, though it has several notable downsides.
First, it is overly complex. Its original formulation is non-convex, and the
optimization requires us to iteratively solve LASSO-equivalent problems until
convergence. The tent-like basis functions add computational cost, and the
number of basis functions used to estimate the coefficient function is
therefore usually restricted (about 230 basis functions were used in previous
works). However, it is often recommended to increase the resolution by
including as many basis functions as computationally feasible, and instead
constrain the model by imposing smoothness (Eilers and Marx, 2021). Moreover,
the finite element framework with tent-like basis functions adds few
advantages. Second, the triangular finite element framework artificially
imposes a constraint on the maximum time lag (i.e., the domain of integration
cannot be negative in Equation 1). In the previous works, each year would be
treated as a separate replicate, independent of all other replicates, which
for our application poses the question of where to start each year, leading
further to problems in determining the cause of streamflow at the beginning of
each replicate where independence between replicates cannot be assumed. For
example, if one were to utilize the previously proposed functional historical
linear models to predict streamflow on day one for the Withlacoochee River in
Figure 2, only one value of rainfall (0mm) could be used, leading to the
impossibility of predicting this year’s peak streamflow. Third, and perhaps
most importantly, as mentioned in future work by Xun, Guan and Cao (2022), the
previous frameworks do not allow for a time dependent time-lag parameter
$\delta(t)$. Other similar methods, with similar issues, can be found in the
distributed lag literature (Almon, 1965; Eisner, 1960; Griliches, 1967;
Trivedi and Lee, 1981; Pesando, 1972; Tanner, 1974; Rushworth et al., 2013).

In our motivating problem, streamflow on some day of the year $t$ is a
function of $\delta(t)$ days of historical rainfall. Both streamflow and
rainfall are measured as daily accumulations in millimeters, so the data is
equally spaced. Ideally, we should be able to learn $\delta(t)$ and the
function $\beta(s,t)$ that different catchments employ to turn rainfall into
streamflow at a daily scale, even in the presence of substantial noise. Given
these circumstances which are not only found in hydrology but also across all
areas of Earth sciences, Eilers (2003) suggests using the identity matrix as a
basis which is equivalent to zero-degree P-splines. This method is referred to
as Whittaker smoothing (Whittaker, 1922), and it vastly simplifies the degree-
three tensor product P-splines used by Rushworth et al. (2013) and the tent
function finite element framework employed by Malfait and Ramsay (2003),
Harezlak et al. (2007), and Xun, Guan and Cao (2022), while not hurting
performance and improving speed (Eilers, 2003; Eilers, Marx and Durbán, 2015;
Eilers and Marx, 2021). This avoids the need to perform costly preprocessing,
removes the need to tune the number of splines (Eilers and Marx, 2021), and
removes the need to assume different functional replicates are independent.
Further, computation is fast because Whittaker smoothing can lead to extremely
sparse matrices (Eilers and Marx, 2021).

This article has three major contributions. First, to the best of our
knowledge, this is the first attempt to estimate a time-dependent time-lag
parameter with historical functional linear models. Second, we simplify and
increase the flexibility of previous historical functional linear model
methods by replacing the tent-like basis functions with the tensor product
Whittaker basis such that higher resolution coefficient and time-lag functions
can be recovered. Finally, we successfully demonstrate that our method can
lead to novel research directions within hydrology by estimating the ground-
truth coefficient function $\beta(s,t)$ and time-dependent time-lag parameter
$\delta(t)$, which can be used in a variety of future applications to
represent the long-term relationship between rainfall and streamflow.

The remainder of this paper is organized as follows. Our method is formally
presented in Section 2 along with its theoretical properties and a sketch of
its algorithm. In Section 3, we conduct several experiments on real streamflow
data and compare the results to expert knowledge. Next, four simulation
experiments in Section 4 assess the quality of parameter estimation at
different levels of noise. Finally, in Section 5, we summarize the work,
explain the limitations of our method and study, and discuss ideas for future
work.

##  2 Method

Let $y_{i}(t)$ be the observed response at time $t\in 1\dots T$ for replicate
$i\in 1...n$, let $x_{i}(t)$ be the corresponding covariate, and let $D$ be
the total number of considered lags (i.e., $\delta(t)\in[0,D-1]$), then the
functional historical linear model with noise term $\epsilon_{i}(t)$ is given
by Equation 2 and visualized in Figure 3.

$y_{i}(t)=\int_{0}^{D-1}\beta(s,t)x_{i}(t-s)\ ds+\epsilon_{i}(t)$ (2)

As with Harezlak et al. (2007) and Xun, Guan and Cao (2022) we assume without
loss of generality that the response and covariate have had their seasonal
signals removed such that both functional variables have mean zero for each
day of the year. Not only does usage of rainfall and streamflow anomalies
instead of the raw signal remove the need for a functional intercept term, it
also removes any potential periodic confounding effects present in both
signals (Moges et al., 2022). In our application, $y_{i}(t)$ is the daily
streamflow during the $i$-th year, and $x_{i}(t)$ is the daily rainfall. We
note that $t-s$ can be negative. When this occurs, we draw from the end of the
previous replicate $i-1$.

Figure 3: The relationships between the response $y$ and the cause $x$ is
shown for D=4.

Let $\phi_{1}(s,t),...,\phi_{K}(s,t)$ denote a sequence of $K$ basis
functions. Given these basis functions and their corresponding basis
coefficients $\mathbb{b}=[b_{1},...,b_{K}]^{T}$, the coefficient function
$\beta(s,t)$ can be denoted as

$\beta(s,t)=\sum_{k=1}^{K}b_{k}\phi_{k}(s,t).$ (3)

Using this basis expansion, we may reformulate Equation 2 as

$\displaystyle y_{i}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{D-1}\sum_{k=1}^{K}b_{k}\phi_{k}(s,t)x_{i}(t-s)\
ds+\epsilon_{i}(t)$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{K}b_{k}\int_{0}^{D-1}\phi_{k}(s,t)x_{i}(t-s)\
ds+\epsilon_{i}(t)$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{K}b_{k}z_{ik}(t)+\epsilon_{i}(t),$

where

$z_{ik}(t)=\int_{0}^{D-1}\phi_{k}(s,t)x_{i}(t-s)\ ds.$

As stated in the introduction, one of our foremost goals is to infer the
dynamic time lag $\delta(t)$ after which the cause has no discernible effect
on the response, i.e.,

$\begin{split}\delta(t)=\max\\{s\\}\hskip 28.45274pt\text{subject to
}\beta(s,t)\neq 0\text{ for a given }t.\end{split}$ (4)

The accuracy in determining $\delta(t)$ is directly related to the resolution
of the basis functions that parameterise $\beta(s,t)$. In each of Harezlak et
al. (2007), Rushworth et al. (2013), and Xun, Guan and Cao (2022), the time
interval is divided into 20 intervals, but when faced with daily data with a
yearly period, this gives a resolution for finding $\delta(t)$ of over two
weeks. Indeed, $\delta(t)$ may change rapidly, but keeping the functional
framework while increasing the number of intervals can be computationally
prohibitive.

Under these circumstances, using the Whittaker basis is the most appropriate
solution (Eilers, 2003; Whittaker, 1922). In this work, we parameterize the
coefficient function $\beta(s,t)$ with tensor products of Whittaker basis
functions $\phi_{k}(s,t)=w_{m}(s)\otimes w_{\ell}(t)$ for
$s\in\\{1,...,D-1\\}$ and $t\in\\{1,...,T\\}$. Here $w_{m}(s)$ and
$w_{\ell}(t)$ denotes the Whittaker basis functions, which are equivalent to
zero-degree B-spline basis functions centered at the data points (Eilers,
2003). If $\phi_{1}(s,t),...,\phi_{K}(s,t)$ are defined as Whittaker basis
functions, the computationally intensive integral term that is found in most
historical functional linear model applications reduces to the data points
themselves.

Suppose the maximum reasonable lagged influence between $x_{i}(t)$ and
$y_{i}(t)$ for all $t$ occurs at lag $s=D-1$, then we must estimate $K=DT$
parameters. In general we expect that $nT=N\ll K$, thus the system is often
underdetermined. This problem cannot be overcome without some assumptions
about the structure of the coefficient function $\beta(s,t)$. Therefore, we
stipulate three major assumptions about the structure of $\beta(s,t)$
throughout this work. First, given that we are at the same lag, coefficients
at consecutive time points should have similar coefficient values (e.g.,
$\beta(s,t)\approx\beta(s,t+1)$). Second, given that we are at the same time
point, coefficients at consecutive lags should have similar values (e.g.,
$\beta(s,t)\approx\beta(s+1,t)$). Finally, we expect $\beta(s,t)$ to be fairly
sparse. Before fitting a model, the extent of the sparsity is unknown, though
we expect the sparsity will be in regions with larger lags. Though not
explicitly stated, these assumptions also persist throughout the previous
works of Harezlak et al. (2007), Xun, Guan and Cao (2022), and Rushworth et
al. (2013). We limit the number of assumptions so that our methods are
flexible and widely applicable (Clark, Kavetski and Fenicia, 2011), though we
could have further assumed that the coefficient function should be non-
negative for our application since additional precipitation cannot reduce
streamflow. Instead of adding this constraint, we keep our model flexible in
this way and consider this as an additional dimension by which we can evaluate
the reliability of our method.

With these three assumptions in mind, we may now impose regularization and
sparsity constraints on $\beta(s,t)$ to reduce the possibility of
overparameterization (Gasparrini et al., 2017). With the first assumption, we
regularize the squared differences of the coefficient function at consecutive
points in time. This can be done with the horizontal first difference penalty
matrix ($D_{H}$), which was first applied to functional historical linear
models in Harezlak et al. (2007) and further used by Xun, Guan and Cao (2022).
We further add periodic boundary conditions to $D_{H}$ in a similar fashion as
Garcia (2010), such that $\beta(s,T)\approx\beta(s,1)$. Likewise, given the
second assumption, we may penalize the squared differences between
coefficients with the same time but with consecutive lag values by introducing
the vertical first difference penalty matrix ($D_{V}$). We additionally
consider a zero top boundary condition on the vertical penalties since lags
past $D-1$ are assumed to have zero influence. The structure of the
coefficients and the penalties are visualized in Figure 4. Further, if we take
$D=2$, $T=3$, and $K=6$, then $D_{V}$, $D_{H}$, and $\beta(s,t)$ are given by:

$\displaystyle D_{V}=\begin{bmatrix}-1&1&0&0&0&0\\\ 0&0&-1&1&0&0\\\
0&0&0&0&-1&1\\\ 0&1&0&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&0&0&0&1\\\
\end{bmatrix}\hskip 28.45274ptD_{H}=\begin{bmatrix}-1&0&1&0&0&0\\\
0&-1&0&1&0&0\\\ 0&0&-1&0&1&0\\\ 0&0&0&-1&0&1\\\ 1&0&0&0&-1&0\\\
0&1&0&0&0&-1\\\ \end{bmatrix}\hskip
28.45274pt\beta(s,t)=\begin{bmatrix}\beta(0,1)\\\ \beta(1,1)\\\ \beta(0,2)\\\
\beta(1,2)\\\ \beta(0,3)\\\ \beta(1,3)\end{bmatrix}.$

Figure 4: The relationships and general structure of the coefficient function
$\beta(s,t)$ are shown for $D=4$ and $T=5$, where $s$ increases along the
vertical axis, and $t$ increases along the horizontal axis. The lines between
nodes indicate that the nodes should have similar values. The partial lines at
the sides of the diagram indicate that there is periodic behavior with a
period of $T+1$ (i.e., $\beta(0,5)$ and $\beta(0,1)$ are connected). The
bottom row of nodes represent the coefficients relating the current time’s
cause to the current time’s response while the next rows represent the
coefficients for consecutive time points in the past.

Most previous works that have developed historical functional linear models
have imposed sparsity in their model (Harezlak et al., 2007; Malfait and
Ramsay, 2003; Xun, Guan and Cao, 2022), however, they all assume that the lag
after which the explanatory feature has no effect on the response is static in
time. This assumption is indeed too strong for many applications, especially
in Earth science. Instead, we impose sparsity with a more flexible approach
via nested group 2-norm thresholding such that a dynamic time lag $\delta(t)$
can be reliably estimated. In statistics, sparsity is usually induced via the
lasso or bridge penalties (Xun, Guan and Cao, 2022; Tibshirani, 1996; Huang,
Horowitz and Ma, 2008). The lasso is known to not have the oracle and sign
consistency properties, meaning it cannot distinguish which coefficients
should be exactly zero according to the data-generating process. This was
shown experimentally in the original lasso paper (Tibshirani, 1996), and it
was shown theoretically in Meinshausen and Yu (2009), Zhao and Yu (2006), Fan
and Li (2001), and Leng, Lin and Wahba (2006). The oracle and consistency
properties of bridge penalties are better than lasso, though they still fail
when the number of parameters exceeds the number of observations and when
there is extensive multicollinearity (Huang, Horowitz and Ma, 2008). Both of
these troublesome properties are expected in problems surrounding our intended
application. As an alternative, Zhao and Yu (2006) suggests using the L0
penalty, though this poses a new problem of non-convexity and NP-hard
computational complexity (Zhao and Yu, 2006; Huo and Ni, 2007). The L0-norm
has the oracle property and is selection consistent, meaning it will select
the correct features which are non-zero in the ground-truth model with a
probability converging to 1 under weaker conditions compared to the lasso
(Zhang and Zhang, 2012; Staerk, Kateri and Ntzoufras, 2018; Kim, Kwon and
Choi, 2012). Due to the strongly hierarchical nature of our model, optimizing
the L0-norm problem changes from $O(2^{K})$ to $O(K)$, thus, for our
application it is computationally tractable and theoretically preferable for
variable selection.

The nested group framework allows us to place each coefficient in a series of
groups according to a hierarchical structure (Zhao, Rocha and Yu, 2009). Let
$A_{0,t},A_{1,t},...,A_{D-1,t}$ be a series of groups defined for each time
point $t\in 1,...,T$. The group $A_{s,t}$ contains the indices of the
coefficients corresponding to $\beta(s,t),...,\beta(D-1,t)$. Notice that the
groups are nested such that $A_{D-1,t}\subset A_{D-2,t}\subset...\subset
A_{0,t}$. The response $Y$ is defined as a column vector of length $N-(D-1)$.
We then define the sparse matrix $\mathbb{Z}\in\mathbb{R}^{N-(D-1)\times K}$
which has $D$ non-zero entries in each row corresponding to the valid entries
of $\beta(s,t)$ given the current time of $y_{i}(t)$. Here, $\mathbb{b}$ is
the vector of $b_{k}\ \forall k\in 1,...,K$. For example, suppose $n=2$,
$D=2$, $T=3$, and $K=6$, then

$\displaystyle\mathbb{Y}=\begin{bmatrix}y_{1}(2)\\\ y_{1}(3)\\\ y_{2}(1)\\\
y_{2}(2)\\\ y_{2}(3)\end{bmatrix}\hskip
14.22636pt\mathbb{Z}=\begin{bmatrix}z_{1,1}(2)&z_{1,2}(2)&z_{1,3}(2)&z_{1,4}(2)&z_{1,5}(2)&z_{1,6}(2)\\\
z_{1,1}(3)&z_{1,2}(3)&z_{1,3}(3)&z_{1,4}(3)&z_{1,5}(3)&z_{1,6}(3)\\\
z_{2,1}(1)&z_{2,2}(1)&z_{2,3}(1)&z_{2,4}(1)&z_{2,5}(1)&z_{2,6}(1)\\\
z_{2,1}(2)&z_{2,2}(2)&z_{2,3}(2)&z_{2,4}(2)&z_{2,5}(2)&z_{2,6}(2)\\\
z_{2,1}(3)&z_{2,2}(3)&z_{2,3}(3)&z_{2,4}(3)&z_{2,5}(3)&z_{2,6}(3)\end{bmatrix}\hskip
14.22636pt\mathbb{b}=\begin{bmatrix}b_{1}\\\ b_{2}\\\ b_{3}\\\ b_{4}\\\
b_{5}\\\ b_{6}\end{bmatrix}.$

Our method first estimates a smooth coefficient function via:

$\min_{\mathbb{b}}||\mathbb{Y}-\mathbb{Z}\mathbb{b}||_{2}^{2}+w_{h}||D_{H}\mathbb{b}||_{2}^{2}+w_{v}||D_{V}\mathbb{b}||_{2}^{2}$
(5)

Then, the group 2-norms are computed and a model is refit using coefficients
with corresponding group norms greater than some threshold $q$ (Equation 6). A
clear and detailed description of our algorithm is provided in Algorithm 1.

$\min_{\mathbb{b}}||\mathbb{Y}-\mathbb{Z}\mathbb{b}||_{2}^{2}+w_{h}||D_{H}\mathbb{b}||_{2}^{2}+w_{v}||D_{V}\mathbb{b}||_{2}^{2}\text{
s.t. }\mathbb{b}_{A_{s,t}}=0\text{ when }||\mathbb{b}_{A_{s,t}}||_{2}^{2}<q$
(6)

0: One response time series $Y=[y_{1},...,y_{N}]$ and one explanatory time
series $X=[x_{1},...,x_{N}]$.

1: Find optimal weights $w_{h},w_{v}$ such that some objective function is
maximized on a validation set.

2: Find optimal $\mathbb{b}$ using the optimal weights $w_{h}$ and $w_{v}$ and
the least squares criterion in Equation 5.

3: Compute the group norms $||\mathbb{b}_{A_{s,t}}||^{2}_{2}$ for all $(s,t)$.

4: Find optimal sparsity threshold $q$.

5: Given $q$, recompute the optimal smoothing weights, $w_{h}$ and $w_{v}$,
and refit the model using all data.

6: return $\beta(s,t)$ and $\delta(t)$

Algorithm 1

##  3 Hydrology Data Application

Since this paper is motivated by problems in hydrology, we use hydrology data
for both our simulated and real data experiments. We select two diverse
catchments to show the versatility of our method. The first catchment we
explore is the Koksilah River which is located in Cowichan, British Columbia,
Canada. The area of this catchment is $236\text{km}^{2}$ with an average
elevation of 461 meters. 11% of precipitation falls as snow and it is an
extremely wet catchment with an aridity index of just 0.37. The second
catchment we explore is the Withlachoochee River which is located in Dade
City, Florida, United States. The area of the catchment is $650\text{km}^{2}$,
and it is low lying with an average elevation of only 42 meters above sea
level. It does not experience any snow and it is a fairly dry catchment with
an aridity index of 1.07. Clearly these two catchments hold significantly
different attributes, making them distinctly interesting case studies
(Janssen, Radić and Ameli, 2021). To obtain rainfall for both these areas, we
use EMDNA (Tang et al., 2021), a high-quality climate dataset that has
complete daily precipitation and temperature values for 1979-2018 at 10km grid
squares across North America. From this data, daily rainfall values are
generated for the Koksilah and Withlacoochee rivers by averaging precipitation
and temperature values across the catchments, then using a temperature
threshold of 0 degrees Celsius, we separate rainfall from snowfall (Figure 5).
Leap days are removed, leaving us with $365*40=14,600$ observations. We
decided that the maximum reasonable lag should be $D=150$ since, no matter the
location, all precipitation in rain-dominated catchments should drain, become
captured by deep groundwater systems, or be evaporated from the catchment
within 150 days (5 months) (Jasechko et al., 2016), leaving a total of
$365*150=54,750$ parameters to be estimated within $\beta(s,t)$. The choice of
maximum lag also alters the total number of available observations since to
predict the response we require the previous 150 days of explanatory data.
Therefore, the final number of observations is 14,451.

Figure 5: The average daily streamflow (red) and rainfall (black) for the
Koksilah River (top) and the Withlacoochee River (bottom) for 1979-2018.

Hyperparameters $w_{h}$ and $w_{v}$ are optimized with Bayesian optimization
given a Gaussian process prior using the GPfit R package (MacDonald, Ranjan
and Chipman, 2015). We start with 30 uniformly random initial points, then
continue by choosing the next point by maximizing the expected increase in
validation set $R^{2}$ for 35 iterations. After experimenting on several
catchments, we noticed that the optimal $w_{h}$ is often much larger than the
optimal $w_{v}$, thus we set our search range as $e^{[10,20]}$ for $w_{h}$ and
$e^{[-5,15]}$ for $w_{v}$.

Choosing the optimal threshold $q$ is often the most difficult step in
Algorithm 1. We do not specifically define how this should be done in the
methods section as how to choose the optimal $q$ is still an open question,
therefore we leave it up to the user and their goals. If the goal is feature
selection, where removing a few important features is not as critical as only
selecting features which are surely important, greater sparsity may be
desired. If instead, a lag should only be removed if it is surely unimportant,
such as in sure independence screening (Fan and Lv, 2008), future users may
desire less sparsity. We aim to satisfy both goals with no specific
preference. In the process of our experiments, we tried several methods.
Naturally, statisticians and data scientists are drawn towards choosing a $q$
that maximizes validation set $R^{2}$. During testing, this method seemed to
give high test set $R^{2}$ scores, however, the final coefficient function had
small values fluctuating around zero and the simulation results were poor,
indicating too little sparsity. The phenomena was also found in Rushworth et
al. (2013) who observed optimizing the AIC leads to too little smoothness,
pushing the authors to use the AIC-optimal smoothness parameter as a lower
bound and artificially increasing this value. The authors believed this point
was important enough to include it as a paramount component of future work.
These results highlight the diverging goals of prediction and inference (Arlot
and Celisse, 2010). A popular alternative method to this is choosing the ”knee
point” in the plot of $q$ versus the whole dataset $R^{2}$ (Satopaa et al.,
2011). This method is popular for unsupervised learning such as clustering
analysis, though for our purposes, it leads to too much sparsity. Instead of
occurring at the middle of the knee, the ground-truth sparsity threshold
seemed to often occur at the beginning of the knee, i.e., ”knee-onset”
(Fermín-Cueto et al., 2020). One promising way to automatically find the
”knee-onset” is maximizing Menger curvature between the first, last, and
current point (Satopaa et al., 2011), though we found that this method was too
sensitive to noise. In all experiments below, we find the ”knee-onset” by
choosing the first of two knees as outlined in Satopaa et al. (2011) and
Fermín-Cueto et al. (2020). All experiments were run in R version 4.2.1.

###  3.1 Analysis of hydrology data

To obtain real response variable data that pairs with our real rainfall data,
daily streamflow data from the Koksilah River was gathered from the HYDAT
database (Albers, 2017), and daily streamflow data from the Withlachoochee
River was taken from CAMELS (Addor et al., 2017). The average daily streamflow
for both rivers is visualized in Figure 5. After performing initial data
exploration and methodological experimentation, we decided to transform
streamflow $Q$ via $log(Q+1)$, since this transform improved predictive
performance and removed patterns and non-Gaussianity from the residual plots.
The estimates of our hyperparameters $w_{h}$ and $w_{v}$ are obtained after
splitting the data with 80% of the data for training and 20% of the data for
validation. The final estimates of $\beta(s,t)$, $\delta(t)$, and the whole
dataset $R^{2}$ are then obtained by training on all available data. We also
evaluate the predictive performance of our model with a test set. To compute
test set $R^{2}$ values, we split the data with 60% for training, 20% for
optimizing the hyperparameters $w_{h}$ and $w_{v}$, and 20% for testing. We
present our results for both rivers below.

The estimated $\hat{\beta}(s,t)$ for the Koksilah River after training on the
entire dataset is shown in Figure 6. The whole dataset $R^{2}$ was 0.74 when
trained on all data. On the test set, the $R^{2}$ was 0.81. In winter,
rainfall frequency and magnitude peaks after a rapid ramp up from September to
October to November (days 260-320), therefore, by December (day 335), the soil
in the catchment is full of water and groundwater is fully connected (Figure
5). When this happens, any additional rain either immediately transmits to
flood waters heading towards the stream, or the additional rainfall quickly
pushes the water already in the soil towards the stream. Further, this
behavior can be seen in Figure 5, where December, January, and February (days
330-70) streamflow correlates strongly with rain falling the previous day.
Therefore, observing a $\delta(t)$ of around 3 days from day 350-60
empirically confirms and supports existing expert knowledge. Two distinct
peaks of maximum lag appear in summer and fall in Figure 6. The peak in summer
at around day 190, with a $\delta(t)$ of about 18 days displays a long term
recession behavior in the catchment. Starting in March, this catchment appears
to enter a drying phase where evapotranspiration and runoff exceed water input
(Figure 5). During the drying phase, rainfall further back in the past begins
to make a larger difference on current streamflow, since the relative amount
of dryness is a significant driver for streamflow production. During August
(around day 230), the soil is completely dry and temperatures are high,
therefore most rainfall contributes to soil wetness and evaporation instead of
streamflow. Starting in mid-September (day 260), the Koksilah River catchment
enters a rapid wetting phase as rainfall quickly becomes more frequent.
Rainfall occurring up to around 40 days in the past begins to have a
significant impact on streamflow, since during this time if the previous days
were dry any additional precipitation will contribute to soil wetness, but if
previous days were wet, additional precipitation will contribute to
streamflow. The wetting phase ends in December (about day 350) when the soil
becomes fully saturated regardless of the between-year variations of
precipitation in the wettest month, November.

Figure 6: The estimated function $\hat{\beta}(s,t)$ for each lag $s$ and time
$t$ for the Koksilah River located in British Columbia from the whole data
period (1979-2018).

For the second catchment, the Withlachoochee River in Florida, the $R^{2}$ was
0.58 on the test dataset. The final whole dataset $R^{2}$ for this catchment
when trained on the entire dataset was 0.54 (Figure 7). It is not surprising
that the whole-dataset $R^{2}$ is often lower than the test set $R^{2}$. This
just indicates that the training set is extremely hard to predict, perhaps due
to high or more frequent extremes which are notoriously hard to predict,
especially with only a linear model. The test set $R^{2}$ results for both
catchments are similar to those seen from LSTMs which are much more complex
and much less interpretable (Hoedt et al., 2021; Anderson and Radić, 2022).

Figure 7: The estimated function $\hat{\beta}(s,t)$ for each lag $s$ and time
$t$ for the Withlachoochee River, located in Florida, after training on the
whole dataset (1979-2018).

Though the Withlachoochee River results presented in Figure 7 look very
different compared to the results from the Koksilah River (Figure 6), the
results can be interpreted in a similar fashion. We first note that this
catchment never reaches a fully saturated phase where rainfall has a huge
impact on streamflow in a short period of time. This is reasonable since this
is a dry catchment with an aridity index of greater than one, meaning that the
potential for evaporation is greater than total precipitation. Starting in May
(days 130-175), the catchment enters a dry phase. By this time of the year,
the catchment has experienced about 6 months of relatively dry weather and
temperatures have begun to rise (Figure 5), therefore in May, any rainfall
that occurs gets evaporated or contributes to soil moisture. By mid to late
June (about day 180), a consistent amount of rainfall usually falls, outpacing
evaporation and streamflow (Figure 5), thus the catchment soils become wet
enough to produce streamflow. Before the catchment can become fully saturated,
in September (day 260), rainfall falls below evaporation and streamflow,
leading to larger lags becoming significant and a drying phase from October to
the next year.

Both real-world studies lead to interesting and meaningful results with good
$R^{2}$ scores on the test set. In both cases, wetting, wet, drying, and dry
periods could be delineated. Further, even though the positivity of the
coefficient function was not imposed as a constraint for our models, in both
cases, the model correctly identified that there is a strong positive
relationship between rainfall and streamflow. A few places with negative
coefficients appeared, but the negative locations were highly limited and the
magnitude of the negative coefficients were extremely small. This is a further
indication that our method is reliable and can correctly identify dominant
processes.

##  4 Simulation Studies

The aim of our simulation study is to quantify and understand our ability to
estimate $\beta(s,t)$ and $\delta(t)$ in a variety of noisy environments.

###  4.1 Evaluation criteria

Both Malfait and Ramsay (2003) and Harezlak et al. (2007) were interested in
model predictions, so in line with their goals, they assessed their methods in
simulation studies by computing $R^{2}$ values between the predicted and
ground-truth response values. In contrast to this, the goal of Xun, Guan and
Cao (2022) was statistical inference, so they compared methods by their
accuracy in estimating $\beta(s,t)$ and $\delta$. More specifically, they
compute the root mean squared error of $\delta$, the percent bias of $\delta$,
and the mean integrated squared error of $\beta(s,t)$. The percent bias of
$\delta$ does not fully capture the accuracy of $\delta$ estimates when
$\delta$ varies over time, and the other metrics are less interpretable
compared to the $R^{2}$ metric used by Malfait and Ramsay (2003) and Harezlak
et al. (2007), thus we introduce three evaluation criteria for our simulation
studies. The first criterion, $\beta$-$R^{2}$, aims to evaluate the accuracy
in estimating the true coefficient function in an interpretable fashion:

$R^{2}(\beta,\hat{\beta})=1-\frac{\int\int\\{\beta(s,t)-\hat{\beta}(s,t)\\}^{2}\
ds\ dt}{\int\int\\{\beta(s,t)-\bar{\beta}\\}^{2}\ ds\ dt}$ (7)

The second criterion, $\delta(t)\text{-bias}$, evaluates the bias appearing in
our estimates of $\delta(t)$. We compute this criterion as the average of
$\hat{\delta}(t)$ across time minus the average of $\delta(t)$ across time.
Our third criterion, $\delta(t)\text{-correlation}$ evaluates how well we can
estimate the inter-period variability in $\delta(t)$ across time. This metric
is simply calculated as the correlation between the ground-truth $\delta(t)$
and the estimated $\hat{\delta}(t)$. Ideally, we would want $\beta$-$R^{2}$
and $\delta(t)\text{-correlation}$ to be close to one, while
$\delta(t)\text{-bias}$ is close to zero.

###  4.2 Simulation scenarios

In this section, we run four simulation studies that aim to quantify our
method’s ability to estimate $\beta(s,t)$ and $\delta(t)$ in noisy
environments. We use the same real rainfall data used in the previous section
to form the covariate matrix $\mathbb{Z}$, then we define four scenarios with
different values for $\beta(s,t)$, $\delta(t)$, and noise levels, such that
known ground-truth values for the response $y_{i}(t)$ can be simulated.
Streamflow is notoriously difficult to predict. Depending on the location of
the catchment and its dynamics, even when using powerful black-box methods
such as long short-term memory neural networks (LSTMs) (Hochreiter and
Schmidhuber, 1997), combined with more predictor variables, the test set
$R^{2}$ can range all the way from zero to one (Ayzel and Heistermann, 2021;
Hoedt et al., 2021; Kratzert et al., 2018). In future scenarios for which our
method may be applied, we hypothesize that inferences from models with
$R^{2}<0.4$ will be difficult. Further, from previous works we know that
interpretable models with $R^{2}>0.8$ will be rare in hydrology. Thus, we
simulate the response vector with additional noise such that if we calculate
the true response $y_{\mathtt{true}}(t)$ from the ground-truth $\beta(s,t)$ we
would produce $R^{2}(y,y_{\mathtt{true}})=0.8$ or
$R^{2}(y,y_{\mathtt{true}})=0.4$, where $R^{2}(y,y_{\mathtt{true}})$ is
defined in Equation 8, $y(t)$ is the simulated data after adding noise to
$y_{\mathtt{true}}(t)$, and $\bar{y}$ is the average of $y(t)$ over $t$.

$R^{2}(y,y_{\mathtt{true}})=1-\frac{\int{\\{y(t)-y_{\mathtt{true}}(t)\\}^{2}\
dt}}{\int{\\{y(t)-\bar{y}\\}^{2}\ dt}}$ (8)

The calculated value of the variance of the added noise for each simulated
scenario is displayed in Table 1.

The four scenarios revolve around the results obtained in the real hydrology
study in Section 3.1. The results from each river were quite different, thus
we take these results from Figures 6 and 7 and set them as our ground-truth
scenarios. Taking these $\beta(s,t)$ values along with our matrix of rainfall
values allows us to compute ground-truth response values. With two different
noise levels, there are four simulation scenarios. For each scenario, we
repeat the simulation 100 times such that stable results and uncertainty
levels can be obtained. Each simulation iteration takes less than 10 minutes
on an Intel i9-9980HK 2.4 GHz processor with 32 GB of RAM, so our methods are
computationally feasible for any modern hardware setup.

###  4.3 Simulation results

Table 1 shows that the introduced method is quite promising and robust. We
consistently observe that regardless of location, the data with less noise
allowed for more accurate and stable inferences of $\beta(s,t)$, however, the
differences are surprisingly small, especially for the Withlacoochee River.
For both rivers, the noise level seemed to make little difference in the bias
of $\delta(t)$ estimates, where both locations saw slight increases in bias
magnitude as we decreased the noise. This suggests that our method for
choosing the threshold $q$ is not optimal. For the Withlachoochee River, the
estimates for $\beta(s,t)$ are consistently accurate with $R^{2}=0.94$ for the
high noise scenario and $R^{2}=0.95$ for the low noise scenario. For both
noise scenarios, the bias was about -6 days while the $\delta(t)$ correlation
was about 0.99. The average results displayed in Table 1 as well as the
individual simulation runs consistently resulted in negative
$\delta(t)$-biases, revealing systematic errors. This indicates that on
average $\delta(t)>\hat{\delta}(t)$, meaning that the estimated coefficient
function is often too sparse and the chosen $q$ should decrease.

Table 1: Summary of the simulation study results. The average (standard deviation) across the 100 replications is shown for each of the evaluation criteria defined in Section 4.1. Location | $R^{2}(y,y_{\mathtt{true}})$ | Noise variance |  $\beta$-$R^{2}$ |  $\delta(t)$-bias |  $\delta(t)$-corr
---|---|---|---|---|---
Koksilah River | 0.4 | 0.244 | 0.940 (0.021) | -0.971 (0.904) | 0.859 (0.094)
Koksilah River | 0.8 | 0.041 | 0.972 (0.006) | -2.16 (0.598) | 0.905 (0.049)
Withlacoochee River | 0.4 | 0.100 | 0.938 (0.010) | -5.559 (0.837) | 0.986 (0.004)
Withlacoochee River | 0.8 | 0.017 | 0.950 (0.005) | -6.478 (0.398) | 0.986 (0.002)

##  5 Discussion and Conclusions

In this work, we build off the iterations of the functional historical linear
model introduced by Malfait and Ramsay (2003), Harezlak et al. (2007),
Rushworth et al. (2013), and Xun, Guan and Cao (2022). We greatly simplify
their formulations while allowing for dynamic sparsity and high resolution
estimates of the coefficient function using the Whittaker basis. In Section 2,
we illustrate our algorithm with several diagrams and introduce the three
minor assumptions we make to reduce the highly overparameterized problem to
one that is tractable. We only assume horizontal smoothness, vertical
smoothness, and that sparsity is concentrated at larger lags to keep our
methods flexible and applicable to a wide variety of applications outside of
our main hydrologically-centered goal. Throughout the work, we are strongly
driven by the goal of accurately estimating rainfall-runoff relationships.
After gathering rainfall and streamflow data from Vancouver Island and
Florida, we estimated the coefficient functions that transfer rainfall to
streamflow at both locations and compared our conclusions with expert
knowledge. We found that the functions have a perfect hydrologically-
significant interpretation where catchments go through four distinct phases:
(1) wetting, (2) wet, (3) drying, (4) dry. Because the Koksilah River
(Vancouver Island) has high seasonal variability in temperature and rainfall,
the four phases can be distinctly parsed. On the other hand, we found that the
Withlachoochee River (Florida) is an extremely dry catchment, therefore it
never enters the wet phase, but all other phases can be inferred. Indeed, none
of these inferences could have been made without the methodological
improvements pertaining to the dynamic sparsity. Finally, in an extensive
simulation study, we found that our simple methodology can successfully
recover the ground-truth coefficient function $\beta(s,t)$ as well as the
dynamic time lag $\delta(t)$ with high accuracy. Though the accuracy for which
we can recover the coefficient function and time-lag clearly depend on the
level of noise, we found that even when the $R^{2}$ is less than or equal to
0.4, accurate conclusions can be made from our inferences.

Several shortcomings of our methods and experiments can be identified. First,
our methods are currently limited to bivariate analysis, though Algorithm 1
could easily be extended to include multiple explanatory features. Second,
though our algorithm can potentially include multiple features, these
additional features may not follow our sparsity and smoothness assumptions.
For example, if we included snowfall as a predictor of streamflow in a snow-
rich catchment we may expect snowfall to take days or weeks or even months to
melt and begin to contribute to streamflow. This would contradict our sparsity
assumption since very low and very high lags would be expected to have zero
coefficients. Third, the way in which we choose sparsity is fairly arbitrary.
Selecting the onset of nonlinearity can be justified by the same theory as the
adjusted-$R^{2}$, since the expected decrease in $R^{2}$ will be linear if we
are removing irrelevant variables, however finding the onset of nonlinearity
in an automatic fashion is still an active area of research. Finally, we
assume throughout the work that the behavior of the catchment does not change
across different years. Although our method does not capture year-to-year
relationship variability, our method can be used to extract multiple
coefficient functions before and after some known change. For example, one
could estimate the rainfall-runoff relationship before and after a large
change in land use or before and after a period of severe climate change.
Focusing now on the limitations of our experiments, we note that our
experiments are only limited to a specific application in hydrology, though we
hypothesize that our methods can work well in other domains. Further, our
simulation study is limited since we provided ground-truths that we know can
be reached from the algorithm. This limitation points back to the second
methodological limitation since our method may only be able to recover certain
types of coefficient functions. If we instead ran a simulation study with a
non-smooth ground-truth, the simulation results would likely be worse.

In future work, one could solve one or more of the above shortcomings. Perhaps
the flexibility of our methods could also be improved by adding adaptive
smoothness in a similar fashion as the adaptive lasso, since different days or
lags could have different smoothness (Zou, 2006; Centofanti et al., 2022; Yang
and Hong, 2017; Yang et al., 2021). Alternatively, extending these methods to
spatial analysis could be interesting. With smooth spatial attribution, one
could understand how impactful a specific grid cell of a cause is on the
response. Further, with sparsity, one could delineate the spatial area that is
impactful on the response versus the areas which have no significant impact on
the response. Our current work takes the perspective of penalized least
squares and sparsity thresholding, but another interesting perspective leading
to different methods could come from a causal inference and conditional
independence framework. For example, finding the correct sparsity using
partial correlation or transfer entropy may be possible (Murari et al., 2018;
Moges et al., 2022), though the high dimensional and low sample size
environment makes this non-trivial. We hope our methods can serve as a
reliable tool for learning from time series data across multiple areas of
science or serve as an inspiration for further methodological development.
This area of research is certainly rich with potential scientific discoveries.

This research was funded by the Collaborative Research Team Project of
Canadian Statistical Sciences Institute (CANSSI) awarded to Ali A. Ameli,
William J. Welch, and Jiguo Cao. Joseph Janssen was supported by Natural
Sciences and Engineering Research Council of Canada (NSERC) PhD Scholarship.
Asad Haris was supported by UBC’s Data Science Institute Postdoctoral funding.
Jiguo Cao’s research was partially supported by the NSERC Discovery grant
(RGPIN-2018-06008). Stefan Schrunner gratefully acknowledges the financial
support from internal funding scheme at Norwegian University of Life Sciences
(project number 1211130114), which financed the international stay at the
University of British Columbia, Canada.

Data and Code. We provide the R code, data, and results for the real data
analysis and simulation studies in the supplementary material.

## References

  * Addor et al. (2017) [author] Addor, NansN., Newman, Andrew JA. J., Mizukami, NaokiN. and Clark, Martyn PM. P. (2017). The CAMELS data set: catchment attributes and meteorology for large-sample studies. Hydrology and Earth System Sciences 21 5293–5313.
  * Albers (2017) [author] Albers, SamS. (2017). tidyhydat: Extract and Tidy Canadian Hydrometric Data. The Journal of Open Source Software 2\. 10.21105/joss.00511
  * Almon (1965) [author] Almon, ShirleyS. (1965). The distributed lag between capital appropriations and expenditures. Econometrica: Journal of the Econometric Society 178–196.
  * Anderson and Radić (2022) [author] Anderson, SamS. and Radić, ValentinaV. (2022). Evaluation and interpretation of convolutional long short-term memory networks for regional hydrological modelling. Hydrology and Earth System Sciences 26 795–825.
  * Arlot and Celisse (2010) [author] Arlot, SylvainS. and Celisse, AlainA. (2010). A survey of cross-validation procedures for model selection.
  * Ayzel and Heistermann (2021) [author] Ayzel, GeorgyG. and Heistermann, MaikM. (2021). The effect of calibration data length on the performance of a conceptual hydrological model versus LSTM and GRU: A case study for six basins from the CAMELS dataset. Computers & Geosciences 149 104708\.
  * Blöschl et al. (2013) [author] Blöschl, GünterG., Bloschl, GunterG., Sivapalan, MurugesuM., Wagener, ThorstenT., Savenije, HubertH. and Viglione, AlbertoA. (2013). Runoff prediction in ungauged basins: synthesis across processes, places and scales. Cambridge University Press.
  * Botter, Bertuzzo and Rinaldo (2010) [author] Botter, GianlucaG., Bertuzzo, EnricoE. and Rinaldo, AndreaA. (2010). Transport in the hydrologic response: Travel time distributions, soil moisture dynamics, and the old water paradox. Water Resources Research 46\.
  * Centofanti et al. (2022) [author] Centofanti, FabioF., Lepore, AntonioA., Menafoglio, AlessandraA., Palumbo, BiagioB. and Vantini, SimoneS. (2022). Adaptive smoothing spline estimator for the function-on-function linear regression model. Computational Statistics 1–26.
  * Chiu and Bittler (1969) [author] Chiu, Chao-LinC.-L. and Bittler, Ronald PR. P. (1969). Linear time-varying model of rainfall-runoff relation. Water Resources Research 5 426–437.
  * Clark, Kavetski and Fenicia (2011) [author] Clark, Martyn PM. P., Kavetski, DmitriD. and Fenicia, FabrizioF. (2011). Pursuing the method of multiple working hypotheses for hydrological modeling. Water Resources Research 47\.
  * Eilers (2003) [author] Eilers, Paul HCP. H. (2003). A perfect smoother. Analytical chemistry 75 3631–3636.
  * Eilers, Marx and Durbán (2015) [author] Eilers, Paul HCP. H., Marx, Brian DB. D. and Durbán, MariaM. (2015). Twenty years of P-splines. SORT: statistics and operations research transactions 39 0149–186.
  * Eilers and Marx (2021) [author] Eilers, Paul HCP. H. and Marx, Brian DB. D. (2021). Practical smoothing: The joys of P-splines. Cambridge University Press.
  * Eisner (1960) [author] Eisner, RobertR. (1960). A distributed lag investment function. Econometrica, Journal of the Econometric Society 1–29.
  * Fan and Li (2001) [author] Fan, JianqingJ. and Li, RunzeR. (2001). Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American statistical Association 96 1348–1360.
  * Fan and Lv (2008) [author] Fan, JianqingJ. and Lv, JinchiJ. (2008). Sure independence screening for ultrahigh dimensional feature space. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70 849–911.
  * Fermín-Cueto et al. (2020) [author] Fermín-Cueto, PaulaP., McTurk, EuanE., Allerhand, MichaelM., Medina-Lopez, EncarniE., Anjos, Miguel FM. F., Sylvester, JoelJ. and Dos Reis, GoncaloG. (2020). Identification and machine learning prediction of knee-point and knee-onset in capacity degradation curves of lithium-ion cells. Energy and AI 1 100006\.
  * Garcia (2010) [author] Garcia, DamienD. (2010). Robust smoothing of gridded data in one and higher dimensions with missing values. Computational statistics & data analysis 54 1167–1178.
  * Gasparrini et al. (2017) [author] Gasparrini, AntonioA., Scheipl, FabianF., Armstrong, BenB. and Kenward, Michael GM. G. (2017). A penalized framework for distributed lag non-linear models. Biometrics 73 938–948.
  * Griliches (1967) [author] Griliches, ZviZ. (1967). Distributed lags: A survey. Econometrica: journal of the Econometric Society 16–49.
  * Harezlak et al. (2007) [author] Harezlak, JaroslawJ., Coull, Brent AB. A., Laird, Nan MN. M., Magari, Shannon RS. R. and Christiani, David CD. C. (2007). Penalized solutions to functional regression problems. Computational statistics & data analysis 51 4911–4925.
  * Hochreiter and Schmidhuber (1997) [author] Hochreiter, SeppS. and Schmidhuber, JürgenJ. (1997). Long short-term memory. Neural computation 9 1735–1780.
  * Hoedt et al. (2021) [author] Hoedt, Pieter-JanP.-J., Kratzert, FrederikF., Klotz, DanielD., Halmich, ChristinaC., Holzleitner, MarkusM., Nearing, Grey SG. S., Hochreiter, SeppS. and Klambauer, GünterG. (2021). Mc-lstm: Mass-conserving lstm. In International Conference on Machine Learning 4275–4286. PMLR.
  * Huang, Horowitz and Ma (2008) [author] Huang, JianJ., Horowitz, Joel LJ. L. and Ma, ShuanggeS. (2008). Asymptotic properties of bridge estimators in sparse high-dimensional regression models. The Annals of Statistics 36 587–613.
  * Huo and Ni (2007) [author] Huo, XiaomingX. and Ni, XueleiX. (2007). When do stepwise algorithms meet subset selection criteria? The Annals of Statistics 870–887.
  * Janssen, Radić and Ameli (2021) [author] Janssen, JosephJ., Radić, ValentinaV. and Ameli, AliA. (2021). Assessment of Future Risks of Seasonal Municipal Water Shortages Across North America. Frontiers in Earth Science 798\.
  * Jasechko et al. (2016) [author] Jasechko, ScottS., Kirchner, James WJ. W., Welker, Jeffrey MJ. M. and McDonnell, Jeffrey JJ. J. (2016). Substantial proportion of global streamflow less than three months old. Nature geoscience 9 126–129.
  * Kim, Şentürk and Li (2011) [author] Kim, KionK., Şentürk, DamlaD. and Li, RunzeR. (2011). Recent history functional linear models for sparse longitudinal data. Journal of statistical planning and inference 141 1554–1566.
  * Kim, Kwon and Choi (2012) [author] Kim, YongdaiY., Kwon, SunghoonS. and Choi, HosikH. (2012). Consistent model selection criteria on high dimensions. The Journal of Machine Learning Research 13 1037–1057.
  * Kratzert et al. (2018) [author] Kratzert, FrederikF., Klotz, DanielD., Brenner, ClaireC., Schulz, KarstenK. and Herrnegger, MathewM. (2018). Rainfall–runoff modelling using Long Short-Term¡? xmltex$\backslash$break?¿ Memory (LSTM) networks. Hydrology and Earth System Sciences 22 6005–6022.
  * Kunz (2020) [author] Kunz, Nadja CN. C. (2020). Towards a broadened view of water security in mining regions. Water Security 11 100079\.
  * Leng, Lin and Wahba (2006) [author] Leng, ChenleiC., Lin, YiY. and Wahba, GraceG. (2006). A note on the lasso and related procedures in model selection. Statistica Sinica 1273–1284.
  * MacDonald, Ranjan and Chipman (2015) [author] MacDonald, BlakeB., Ranjan, PritamP. and Chipman, HughH. (2015). GPfit: An R package for fitting a Gaussian process model to deterministic simulator outputs. Journal of Statistical Software 64 1–23.
  * Malfait and Ramsay (2003) [author] Malfait, NicoleN. and Ramsay, James OJ. O. (2003). The historical functional linear model. Canadian Journal of Statistics 31 115–128.
  * Meinshausen and Yu (2009) [author] Meinshausen, NicolaiN. and Yu, BinB. (2009). Lasso-type recovery of sparse representations for high-dimensional data. The annals of statistics 37 246–270.
  * Moges et al. (2022) [author] Moges, EdomE., Ruddell, Benjamin LB. L., Zhang, LiangL., Driscoll, Jessica MJ. M. and Larsen, Laurel GL. G. (2022). Strength and Memory of Precipitation’s Control Over Streamflow Across the Conterminous United States. Water Resources Research 58 e2021WR030186.
  * Murari et al. (2018) [author] Murari, AndreaA., Lungaroni, MicheleM., Peluso, EmmanueleE., Gaudio, PasqualeP., Lerche, ErnestoE., Garzotti, LucaL., Gelfusa, MichelaM. and Contributors, JETJ. (2018). On the use of transfer entropy to investigate the time horizon of causal influences between signals. Entropy 20 627\.
  * Pesando (1972) [author] Pesando, James EJ. E. (1972). Seasonal variability in distributed lag models. Journal of the American Statistical Association 67 311–312.
  * Rushworth et al. (2013) [author] Rushworth, Alastair MA. M., Bowman, Adrian WA. W., Brewer, Mark JM. J. and Langan, Simon JS. J. (2013). Distributed lag models for hydrological data. Biometrics 69 537–544.
  * Satopaa et al. (2011) [author] Satopaa, VilleV., Albrecht, JeannieJ., Irwin, DavidD. and Raghavan, BarathB. (2011). Finding a” kneedle” in a haystack: Detecting knee points in system behavior. In 2011 31st international conference on distributed computing systems workshops 166–171. IEEE.
  * Schuite et al. (2019) [author] Schuite, JonathanJ., Flipo, NicolasN., Massei, NicolasN., Rivière, AgnèsA. and Baratelli, FulviaF. (2019). Improving the spectral analysis of hydrological signals to efficiently constrain watershed properties. Water Resources Research 55 4043–4065.
  * Staerk, Kateri and Ntzoufras (2018) [author] Staerk, ChristianC., Kateri, MM. and Ntzoufras, II. (2018). Adaptive subspace methods for high-dimensional variable selection, PhD thesis, Universitätsbibliothek der RWTH Aachen.
  * Tang et al. (2021) [author] Tang, GuoqiangG., Clark, Martyn PM. P., Papalexiou, Simon MichaelS. M., Newman, Andrew JA. J., Wood, Andrew WA. W., Brunet, DominiqueD. and Whitfield, Paul HP. H. (2021). EMDNA: An Ensemble Meteorological Dataset for North America. Earth System Science Data 13 3337–3362.
  * Tanner (1974) [author] Tanner, J ErnestJ. E. (1974). Variable distributed lags and forecasting non-residential construction. Canadian Journal of Economics 642–654.
  * Tibshirani (1996) [author] Tibshirani, RobertR. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological) 58 267–288.
  * Trivedi and Lee (1981) [author] Trivedi, Pravin KP. K. and Lee, BMSB. (1981). Seasonal variability in a distributed lag model. The Review of Economic Studies 48 497–505.
  * Vink and Robbins (2012) [author] Vink, SueS. and Robbins, SarahS. (2012). Mine water management in variable climate regimes. In Presented at the International Mine Water Association Symposium. pp. 20A–20I.
  * Whittaker (1922) [author] Whittaker, Edmund TE. T. (1922). On a new method of graduation. Proceedings of the Edinburgh Mathematical Society 41 63–75.
  * Xun, Guan and Cao (2022) [author] Xun, XiaoleiX., Guan, TianyuT. and Cao, JiguoJ. (2022). Sparse estimation of historical functional linear models with a nested group bridge approach. Canadian Journal of Statistics. https://doi.org/10.1002/cjs.11747
  * Yang and Hong (2017) [author] Yang, LianqiangL. and Hong, YongmiaoY. (2017). Adaptive penalized splines for data smoothing. Computational Statistics & Data Analysis 108 70–83.
  * Yang et al. (2021) [author] Yang, LianqiangL., Ding, MengzhenM., Hong, YongmiaoY. and Wang, XuejunX. (2021). Estimating functions and derivatives via adaptive penalized splines. Communications in Statistics-Simulation and Computation 50 2054–2071.
  * Zhang and Zhang (2012) [author] Zhang, Cun-HuiC.-H. and Zhang, TongT. (2012). A general theory of concave regularization for high-dimensional sparse estimation problems. Statistical Science 27 576–593.
  * Zhao, Rocha and Yu (2009) [author] Zhao, PengP., Rocha, GuilhermeG. and Yu, BinB. (2009). The composite absolute penalties family for grouped and hierarchical variable selection. The Annals of Statistics 37 3468–3497.
  * Zhao and Yu (2006) [author] Zhao, PengP. and Yu, BinB. (2006). On model selection consistency of Lasso. The Journal of Machine Learning Research 7 2541–2563.
  * Zou (2006) [author] Zou, HuiH. (2006). The adaptive lasso and its oracle properties. Journal of the American statistical association 101 1418–1429.