*2.1. Mathematical Model*

Let {**<sup>X</sup>***i*,**Y***i*}*ni*=<sup>1</sup> be a set of observations of a stochastic process, **X** = -*<sup>X</sup>*<sup>1</sup>(*t*),..., *Xp*(*t*), where *Xj*(*t*) ∈ *L*2 [0, *T*] , *j* = 1, ... , *p*, are predictor covariates and **Y** = (*<sup>Y</sup>*1,*Y*2), with *Yj* ∈ R, a response variable. In this context, the following bivariate location-scale model [40,41] is assumed

$$
\pi \left( \begin{array}{c} Y\_1 \\ Y\_2 \end{array} \right) = \left( \begin{array}{c} \mu\_1(\mathbf{X}) \\ \mu\_2(\mathbf{X}) \end{array} \right) + \mathbf{\mathcal{L}}^{1/2}(\mathbf{X}) \left( \begin{array}{c} \varepsilon\_1 \\ \varepsilon\_2 \end{array} \right) \tag{1}
$$

where **Σ**1/2(**X**) represents the Cholesky decomposition of the variance-covariance matrix **Σ**(**X**)

$$\Sigma(\mathbf{X}) = \begin{pmatrix} \sigma\_1^2(\mathbf{X}) & \sigma\_{12}(\mathbf{X}) \\ \sigma\_{12}(\mathbf{X}) & \sigma\_2^2(\mathbf{X}) \end{pmatrix} \tag{2}$$

so that Var(**Y**|**X**) = **Σ**(**X**) = **Σ1**/**<sup>2</sup>**(**X**) **Σ1**/**<sup>2</sup>**(**X**)**<sup>T</sup>**. To guarantee the model identification in (1), the bivariate residuals (*<sup>ε</sup>*1,*ε*2) are assumed to be independent of the covariates, with zero mean, unit variance, and zero correlation. Despite we do not assume any distribution for the error term, within the framework of functional data analysis this work might be addressed under the assumption of other structures for error distribution: generalized Gauss-Laplace distribution that relax the constrictive assumption of the normal distribution errors [42], generalized linear mixed models (GLMMs) [38] to estimate random effects and dependent (temporal or spatial) errors, and generalized additive models for location, scale and shape (GAMLSS) [43] to model the dynamically variable distribution, considering skewness and kurtosis.

We define the unconditionally probabilistic region for the errors (*<sup>ε</sup>*1,*ε*2) as

$$\varepsilon\_{\mathsf{T}}(k) = \{ (\varepsilon\_1, \varepsilon\_2) \in \mathbb{R}^2 | f(\varepsilon\_1, \varepsilon\_2) \ge k \}$$

*f* being the density function of the bivariate residuals (*<sup>ε</sup>*1,*ε*2) and *k* the *<sup>τ</sup>*−quantile of *f*(*<sup>ε</sup>*1,*ε*2). Then, for a given **X**, we define a conditional *<sup>τ</sup>th*- uncertainty region for (*<sup>Y</sup>*1,*Y*2) containing *τ*% of the observations as:

$$R\_{\mathsf{T}}(\mathsf{X}) = \begin{pmatrix} \mu\_1(\mathsf{X}) \\ \mu\_2(\mathsf{X}) \end{pmatrix} + \mathsf{L}^{1/2}(\mathsf{X})\varepsilon\_{\mathsf{T}} \tag{3}$$
