Varying-Coefficient Additive Models with Density Responses and Functional Auto-Regressive Error Process

1. Introduction

Density, or more generally, distributional data appear increasingly often in various real-world research domains. Specific examples include the distributions of cross-sectional or intraday stock returns [9,15], mortality densities [13], and distributions of intrahub connectivity in neuroimaging [11,14]. In many applications, the density curves are observed consecutively over time, which we call in the article a density time series. A motivating example is shown in Figure 1, (a) shows density time series of global mortality rate (‰) over an interval of 100 days from January 22, 2020 to April 15, 2021, while Figure (b) is an alternative look with plotting densities at three different days. In this article, we consider a regression model where density time series, served as response, is coupled with scalar predictors.

Unlike conventional functional data, density data do not form a linear space due to inherent constraints, such as nonnegativity and the requirement to integrate to one. These characteristics present significant challenges when attempting to directly apply functional data analysis techniques to random densities. Several studies have addressed this issue and proposed various solution methods, which can be broadly categorized into two types. One natural approach to overcoming the nonlinear limitations of density functions is to map them into a Hilbert space using a suitable transformation method. In this direction, [12] proposed two continuous and invertible transformations, which are log hazard transformation and log quantile density (LQD) transformation, to map probability densities to an unrestricted space of square integrable functions. [8] map densities by LQD transformation as in [12] to unrestricted square integrable functions, model and fit the responses with an additive functional-to-scalar regression. In order to forecast density functions derived from cross-sectional and intraday returns, [9] proposed two approaches based on compositional data analysis and modified log quantile transformation, combining the FPC representation and exponential smoothing forecasting method. The second way involves choosing appropriate metrics and adopting a geometric approach. For instance, [17] adopted the infinite-dimension version of the Aitchison geometry to construct a density-on-scalar linear regression model in Bayes–Hilbert spaces, while [13] discussed Fréchet regression in a general metric space with the Wasserstein metric. [5] utilized the geometry of tangent bundles of the Wasserstein space with the Wasserstein metric to propose distribution-on-distribution regression. They also studied an extension to autoregressive models of order one for distribution-valued time series. With different methodology, [19] also discussed order-p Wasserstein autoregressive model.

Denote $\mathcal{F}$ as the space of density functions f defined on a common support $\mathcal{U}$. Without of generality, we assume that $\mathcal{U}=[0,1]$. Given the transformation $\Psi :\mathcal{F}\to L2$, for any given $\mathbf{X}\in {\mathbb{R}}^d$, the conditional Fréchet mean of the random density f is defined as

$$\mu (·|\mathbf{X})=arg\underset{d\in \mathcal{F}}{min}E\left(\left|\right|\Psi \left(f\right)-{\Psi \left(d\right)\left|\right|}_2^2\right|\mathbf{X}),$$

where the expectation E represents the joint distribution of $(\mathbf{X},f)$.

This is equivalent to the following equation

$$\Psi \left(\mu (·|\mathbf{X})\right)\left(u\right)=E\left(\Psi \left(f\right)\left(u\right)|\mathbf{X}\right),0\le u\le 1,$$

leading to the fact that

$$\mu \left(s\right|\mathbf{X})={\Psi }^{-1}\left(E(\Psi (f\left)\right(u\left)\right|\mathbf{X})\right)\left(s\right),0\le s\le 1.$$

Data that we consider in this article are the density time series $d_t$ coupled with scalar predictors $({\mathbf{X}}_t,{\mathbf{Z}}_t)$. To deal with the density functions, we adopt the LQD transformation $\Psi :\mathcal{F}\to L_2$, where $\mathcal{F}$ is the class of density function d satisfying ${\int }_R{u}^{2}d\left(u\right)du<\infty$. For each $d_t\in \mathcal{F}$, let $Q_t\left(u\right)$ be the quantile function corresponding to the cumulative distribution function $F_t\left(y\right)$ of $d_t\left(y\right)$ with support on [0,1], and $q_t\left(u\right)$ be the quantile density function, i.e., $q_t\left(u\right)=Q_t^{\prime }\left(u\right)={\displaystyle \frac{d}{du}}{F}_t^{-1}\left(u\right)$ for $u\in [0,1]$, then the LQD transformation of $d_t$ is defined as

$$\Psi \left(d_t\right)\left(u\right)=log\left({\displaystyle \frac{d}{du}}{F}_t^{-1}\left(u\right)\right),u\in [0,1].$$

In this article, a varying-coefficient additive model with functional auto-regressive error process is proposed to estimate $E(\Psi \left(d_t\right)|{\mathbf{X}}_t,{\mathbf{Z}}_t)$. Then $d_t$ can be written as

$$d_t={\Psi }^{-1}\left(E(\Psi \left(d_t\right)|{\mathbf{X}}_t,{\mathbf{Z}}_t)\right)+{\delta }_{t1},$$

where ${\delta }_{t1}$ is the error from regression.

Except for ${\delta }_{t1}$, there is usually another source of error, which comes from the estimation of $d_t$, i.e., ${\widehat{d}}_t=d_t+{\delta }_{t2}$. That is because in most practical applications, instead of the density $d_t$ itself, only $n_t$ random samples generated from each $d_t$ can be observed. Namely, $Y_{t1},\dots ,Y_{t{n}_t}\sim d_t$. In this article, we assume that $n_t=n$. Following [12], $d_t$ is estimated by the modified kernel smoothing method, i.e,

$${\widehat{d}}_t\left(y\right)=\sum _{i=1}^n\mathcal{K}\left({\displaystyle \frac{y-Y_{ti}}{h}}\right)w(y,h)/\sum _{i=1}^n{\int }_0^1\mathcal{K}\left({\displaystyle \frac{s-Y_{ti}}{h}}\right)w(s,h)ds,$$

where the weight function $w(y,h)$ is defined so as to rectify the boundary effects as follows:

$$w(y,h)={\left({\int }_{-y/h}^1\mathcal{K}\left(u\right)du\right)}^{-1}{I}_{y\in [0,h)}+{\left({\int }_{-1}^{(1-y)/h}\mathcal{K}\left(u\right)du\right)}^{-1}{I}_{y\in (1-h,1]}+I_{y\in [h,1-h]},$$

in which $\mathcal{K}$ is the kernel with bandwidth $h<1/2$ and bounded variation that is symmetric about 0 such that ${\int }_0^1\mathcal{K}\left(u\right)du>0$, ${\int }_{\mathcal{R}}\left|u\right|\mathcal{K}\left(u\right)du$, ${\int }_{\mathcal{R}}{\mathcal{K}}^2\left(u\right)du$, and ${\int }_{\mathcal{R}}\left|u\right|{\mathcal{K}}^2\left(u\right)du$ are finite. Thus, when we fit the regression model with ${\widehat{d}}_t$ as the observation of $d_t$, we have

$${\widehat{d}}_t={\Psi }^{-1}\left(E(\Psi \left(d_t\right)|{\mathbf{X}}_t,{\mathbf{Z}}_t)\right)+{\delta }_{t1}+{\delta }_{t2}.$$

The contribution of this article lies in integrating density time series modeling with a functional autoregressive error process. A common assumption in regression modeling is the independence of random errors. However, data collected over time often exhibit serial dependence, rendering this assumption potentially inappropriate. By incorporating autoregressive structures into the error process, this work accounts for the temporal correlation inherent in density time series, thereby improving the model’s flexibility and accuracy. For example, daily price curves or liquidity demand-supply curves in economic and financial applications or climatological curves describing the natural phenomenon [1,3,4,6] all change over time. In such cases, the daily fluctuations in the collected data are heavily influenced by previous values, potentially resulting in what is known as sequence dependence. Numerous empirical studies have demonstrated that neglecting the temporal dependencies between observations can lead to biased statistical results and flawed inferences. Consequently, it is crucial to incorporate autoregressive error structures that explicitly account for the autocorrelation present within the data. This approach ensures a more accurate representation of the underlying temporal dynamics in regression modeling. While functional autoregressive (FAR) models have been extensively studied in the literature [1,2,3,4,6,10,18], there is a notable gap in the literature regarding the use of FAR processes to model the random error component in functional regression models. In this article, to improve the accuracy of density time series predictions, we integrate the FAR error structure within a varying coefficient additive regression model. This novel approach allows for better handling of temporal dependencies in the error process, ultimately enhancing the model’s predictive performance.

The remainder of this article is organized as follows. Section 2 outlines the methodology for constructing a varying coefficient additive model with a density response, incorporating a functional autoregressive (FAR) error process. In Section 3, we introduce a three-step procedure for estimating the bivariate varying-coefficient components within the model. Section 4 provides the theoretical results and corresponding inferences derived from the model. Monte Carlo simulations are conducted in Section 5 to demonstrate the efficiency and robustness of the proposed approach. In Section 6, we showcase the practical applicability of the model through analyses of COVID-19 data and U.S. income data. Finally, a discussion of the findings is provided in Section 7, with detailed proofs of the theoretical results included in the Appendix.

2. Model Setup

In this article, we focus on the modeling of density response. Due to the constraints of density function which is nonnegative and integrated to 1, we consider the functions after LQD transformation.

In the article, the main goal is to estimate $E(\Psi \left(d_t\right)\left(u\right)|\mathbf{x})$ based on the transfer of density function as

$$E(\Psi \left(d_t\right)\left(u\right)|\mathbf{x})=\sum _{m=1}^k{z}_{t,m}{g}_m(u,x_{t,m}),0\le u\le 1,$$

(1)

leading to our proposed varying-coefficient additive models with density responses and functional auto-regressive $FAR\left(p\right)$ error process (DVCA-FAR):

$$f_t\left(u\right)=\Psi \left(d_t\right)\left(u\right)=\sum _{m=1}^k{z}_{t,m}{g}_m(u,x_{t,m})+{\epsilon }_t\left(u\right),0\le u\le 1,1\le t\le T,$$

(2)

where ${\epsilon }_t\left(u\right)$ is the $FAR\left(p\right)$ process

$${\epsilon }_t\left(u\right)=\int {\gamma }_1(s,u){\epsilon }_{t-1}\left(s\right)ds+\dots +\int {\gamma }_p(s,u){\epsilon }_{t-p}\left(s\right)ds+e_t\left(u\right).$$

(3)

In the above model, the random density $d_t(·)\in \mathcal{F}$ serves as a response, and $\Psi :\mathcal{F}⟶L_2$ is the LQD transformation. Each density is coupled with the $k-$dimensional covariates ${\mathbf{x}}_t={(x_{t,1},\dots ,x_{t,k})}^{\tau }$ and ${\mathbf{z}}_t={(z_{t,1},\dots ,z_{t,k})}^{\tau }$ with supports $S_x$ and $S_z$ respectively. Without loss of generality, we take $S_x=S_z=[0,1]$. In this article, covariate ${\mathbf{x}}_t$ can be ${\mathbf{z}}_t$ or $t/T$.

The bivariate functions $g_m(·,x_m)$ quantify the effect of $\mathbf{z}$, ${\gamma }_l(·,·)$ are smooth functions satisfying $\int \int {\gamma }_l^2(s,u)duds<\infty$. Moreover, the error $e_t\left(u\right)$ is an independent and identically distributed random process satisfying $E\left(e_t\left(u\right)\right)=0$ and $Cov(e_t\left(u\right),e_t\left(s\right))={\sigma }_t^2(u,s)$.

With the estimated density $\widehat{d}$, let ${\widehat{f}}_t=\Psi \left({\widehat{d}}_t\right)$. Then DVCA-FAR model can be written as

$${\widehat{f}}_t\left(u\right)=\sum _{m=1}^k{z}_{t,m}{g}_m(u,x_{t,m})+{\epsilon }_t\left(u\right)+{\epsilon }_{f_t}\left(u\right),0\le u\le 1,$$

(4)

where ${\epsilon }_{f_t}\left(u\right)$ is random error generated from the transformation of estimated density.

3. Three-Step Estimation Methodology

We propose a three-step estimation procedure. First, the B-spline smoothing method is employed to obtain the initial estimate of the bivariate varying coefficients, disregarding the error structure. Second, using the initial estimate and the response, we estimate the error process. The functional autoregressive (FAR) process is then estimated using the sequential test proposed by [10]. Finally, after removing the estimated FAR error, the optimal estimate is derived using the spline method.

3.1. Initial Estimation of Bivariate Varying-Coefficient Function

First, without considering error structure, we adopt the spline method to obtain the initial estimations of $g_m(u,x_m)$, $m=1,\dots ,k,$ regardless of the error process, by the spline series approximation method.

Let $\{B_0\left(u\right),\dots ,B_{N_0}\left(u\right)\}$ be the set of B-spline basis functions for u of order q with $L_0$ interior knots, where $N_0+1=L_0+q$. Similarly, for each $m=1,\dots ,k$, let $\{B_{0,m}\left(x_m\right),\dots ,B_{N_m,m}\left(x_m\right)\}$ be the set of B-spline basis functions of order q for $x_m$ ($m=1,\dots ,k$) with $L_m$ interior knots, where $N_m+1=L_m+q$. Denote $b_{j,m}^{*}\left(x_m\right)$ as the normalization of $B_{j,m}\left(x_m\right)$ and $b_r\left(u\right)=N_0^{1/2}{B}_r\left(u\right)$ as the scaled version of $B_r\left(u\right)$.

Using these basis functions, we define the tensor product of the B-spline basis as

$$b_{r,j,m}(u,x_m)=b_r\left(u\right)b_{j,m}^{*}\left(x_m\right),1\le r\le N_0,1\le j\le N_m,1\le m\le k.$$

The spline approximation of $g_m(u,x_m)$ is given by

$$g_m(u,x_m)\approx \sum _{r=1}^{N_0}\sum _{j=1}^{N_m}{\lambda }_{r,j,m}{b}_{r,j,m}(u,x_m),1\le m\le k.$$

With the least square method, the estimation of $g_m(u,x_m)$ is

$${\tilde{g}}_m(u,x_m)=\sum _{r=1}^{N_0}\sum _{j=1}^{N_m}{\tilde{\lambda }}_{r,j,m}{b}_{r,j,m}(u,x_m),1\le m\le k,$$

(5)

where $\tilde{\lambda }={({\tilde{\lambda }}_{1,1,1},\dots ,{\tilde{\lambda }}_{N_0,N_k,k})}^{\tau }$ is a $\left(N_0{\sum }_{m=1}^k{N}_m\right)$-dimensional vector satisfying

$$\tilde{\lambda }=arg\underset{\lambda }{min}\sum _{t=1}^T\sum _{i=1}^n{\left[{\widehat{f}}_t\left(u_i\right)-\sum _{m=1}^k{z}_{t,m}\sum _{r=1}^{N_0}\sum _{j=1}^{N_m}{\lambda }_{r,j,m}{b}_{r,j,m}(u_i,x_{t,m})\right]}^2.$$

(6)

Theorem 1 shows that the initial estimation ${\tilde{g}}_m(u,x_m)$ is uniformly consistent.

3.2. Estimation of FAR Error Process

With the initial estimation ${\tilde{g}}_m(u,x_m)$, we now estimate the FAR error process. To this end, let

$${\tilde{\epsilon }}_t\left(u\right)=f_t\left(u\right)-\sum _{m=1}^k{z}_{t,m}{\tilde{g}}_m(u,x_{t,m}),1\le t\le T.$$

(7)

Denote the additive function ${\rho }_t\left(u\right)={\sum }_{l=1}^p\int {\gamma }_l(s,u){\epsilon }_{t-l}\left(s\right)ds$, thus the $FAR\left(p\right)$ error process (3) can be written as ${\epsilon }_t\left(u\right)={\rho }_t\left(u\right)+e_t\left(u\right)$.

Let $\{B_0\left(u\right),B_2\left(u\right),\dots ,B_N\left(u\right)\}$ be the set of B-spline basis functions of order q with L interior knots, where $N+1=L+q$. Define the tensor product of the B-spline basis as

$$b_{r,j}(u,s)=B_r\left(u\right)B_j\left(s\right),1\le r,j\le N.$$

Then, the spline approximations of the functions ${\gamma }_l(·,·)$ is given by

$${\gamma }_l(s,u)=\sum _{r=1}^N\sum _{j=1}^N{\mu }_{r,j,l}{b}_{r,j}(u,s),1\le l\le p.$$

The estimation of the $p{N}^2$-dimensional vector $\mu ={({\mu }_{1,1,1},\dots ,{\mu }_{N,N,p})}^{\tau }$ can be obtained by minimizing the following quadratic loss function:

$$\widehat{\mu }=arg\underset{\mu }{min}\sum _{t=p+1}^T\sum _{i=1}^n{\left[{\tilde{\epsilon }}_t\left(u_i\right)-\sum _{l=1}^p\sum _{r=1}^N\sum _{j=1}^N{\mu }_{r,j,l}\int b_{r,j}(u_i,s){\tilde{\epsilon }}_{t-l}\left(s\right)ds\right]}^2.$$

Then, the corresponding estimation of ${\gamma }_l(·,·)$ and the additive function ${\rho }_t\left(u\right)$ are given by

$${\widehat{\gamma }}_l(s,u)=\sum _{r=1}^N\sum _{j=1}^N{\widehat{\mu }}_{r,j,l}{b}_{r,j}(u,s),1\le l\le p,$$

and

$${\widehat{\rho }}_t\left(u\right)=\sum _{l=1}^p\sum _{r=1}^N\sum _{j=1}^N{\widehat{\mu }}_{r,j,l}\int b_{r,j}(u,s){\widehat{\epsilon }}_{t-l}\left(s\right)ds,0\le u\le 1,p+1\le t\le T,$$

(8)

respectively.

In practice, the order p of error process is unknown. To solve this problem, we utilize a sequential test by [10] to identify the order p. The detailed identification procedure is given in Section 3.4.2.

3.3. Improved Estimation of Bivariate Varying-Coefficient Function

By removing the estimated error process (8) from the response function $f_t\left(u\right)$ and repeating the same procedure as in Section 3.1, we obtain an improved estimation of $g_m(u,x_m)$.

To do so, we first denote

$$f_t^c\left(u\right)=f_t\left(u\right)-\sum _{l=1}^p\int {\gamma }_l(s,u){\epsilon }_{t-l}\left(s\right)ds,0\le u\le 1,p+1\le t\le T,$$

and the estimation of $f_t^c\left(u\right)$ as

$${\widehat{f}}_t^c\left(u\right)=f_t\left(u\right)-\sum _{l=1}^p\int {\widehat{\gamma }}_l(s,u){\widehat{\epsilon }}_{t-l}\left(s\right)ds,0\le u\le 1,p+1\le t\le T.$$

On the other side,

$$f_t^c\left(u\right)=\sum _{m=1}^k{z}_{t,m}{g}_m(u,x_{t,m})+e_t\left(u\right),0\le u\le 1.$$

Therefore, following the same procedure as the first step, the improved spline approximation estimations are given by

$${\widehat{g}}_m(u,x_m)=\sum _{r=1}^{N_0}\sum _{j=1}^{N_m}{\widehat{\lambda }}_{r,j,m}{b}_{r,j,m}(u,x_m),1\le m\le k,$$

(9)

where $\widehat{\lambda }={({\widehat{\lambda }}_{1,1,1},\dots ,{\widehat{\lambda }}_{N_0,N_k,k})}^T$ is a $\left(N_0{\sum }_{m=1}^k{N}_m\right)$-dimensional vector satisfying

$$\widehat{\lambda }=arg\underset{\lambda }{min}\sum _{t=1}^T\sum _{i=1}^n{\left[{\widehat{f}}_t^c\left(u_i\right)-\sum _{m=1}^k{z}_{t,m}\sum _{r=1}^{N_0}\sum _{j=1}^{N_m}{\lambda }_{r,j,m}{b}_{r,j,m}(u_i,x_{t,m})\right]}^2.$$

(10)

Theorem 2 and 3 show the uniform convergence and asymptotic normality of the improved estimation ${\widehat{g}}_m(u,x_m)$. Besides, simulation studies in Section 5 show that the improved estimation is more efficient than the initial estimation ${\tilde{g}}_m(u,x_m)$.

3.4. Implementation

3.4.1. Selection of Bandwidth

In empirical applications, prior to modeling, it is necessary to estimate the density, which involves selecting the appropriate bandwidth for the modified kernel estimation. In this section, we employ the leave-one-out cross-validation method to determine the optimal bandwidth. Specifically, the bandwidth h is chosen by minimizing the following mean squared error (MSE):

$$CV\left(h\right)={\displaystyle \frac{1}{nT}}\sum _{t=1}^T\sum _{i=1}^n{[d_t\left(y_{ti}\right)-{\widehat{d}}_t\left(y_{ti}\right)]}^2,$$

where for each $i=1,\dots ,n$, ${\widehat{d}}_t\left(y_{ti}\right)$ is the estimate of $d_t\left(y_{ti}\right)$ with bandwidth h obtained by using observations from the t-th section other than the i-th one.

3.4.2. Identifying the Order of the FAR Process

In this section, we utilize the order determination procedure proposed by [10] to identify th order p of FAR error process. The main idea is to represent the FAR process as a fully functional linear model with dependent regressors and construct a sequential test procedure.

Given a sequence of hypotheses:

$$H_{0,p}:\left\{{\epsilon }_t\right\}\mathrm{are}FAR\left(p\right)vs{H}_{a,p+1}:\left\{{\epsilon }_t\right\}\mathrm{are}FAR(p+1),p=0,1,2,\dots ,$$

Here, $FAR\left(0\right)$ means independent and identically distributed process.The sequential test begins with $p=0$ and terminated if $H_{0,p}$ is not rejected, the order of the process is then identified as p. See [10] for more relevant details and further explanation.

To construct the test statistic, define ${\eta }_j\left(s\right)={\sum }_{l=1}^p{\tilde{\epsilon }}_{j-l}(sp-(l-1))I_l\left(s\right),\phi (s,u)=p{\sum }_{l=1}^p{\gamma }_l(sp-(l-1),u)I_l\left(s\right),$ where $I_l$ as the indicator function of the interval $\left[\right(l-1)/p,l/p]$. Denote the orthonormal basis of $L^2$ as $\{{\widehat{x}}_j,1\le j\le T\}$ from the eigenfunctions of ${\widehat{C}}_{\eta }(s,u)={\displaystyle \frac{1}{T}}{\sum }_{j=1}^T({\eta }_j\left(s\right)-\overline{\eta }\left(s\right))({\eta }_j\left(u\right)-\overline{\eta }\left(u\right)),$ with the corresponding decreasingly ordered eigenvalues ${\widehat{\lambda }}_j$, where $\overline{\eta }$ is the mean function of ${\eta }_j$. In the proposed method, only the first $q_{\eta }$ eigenfunction/eigenvalue pairs are used. Moreover, define $\{{\widehat{y}}_j,1\le j\le T\}$ and $q_{\pi }$ analogously to $\left\{{\widehat{x}}_j\right\}$ and $q_{\eta }$ for the response functions denoted as ${\pi }_j$.

For the product space $L^2([0,1]\times [0,1])$, denote $\eta (j,k)=<{\eta }_j,{\widehat{x}}_k>$, $\pi (j,m)=<{\pi }_j,{\widehat{y}}_m>$, $\psi (k,m)=<\phi ,{\widehat{x}}_k\otimes {\widehat{y}}_m>$. Let $\pi ={\left[\pi (j,m)\right]}_{T\times q_{\pi }},\eta ={\left[\eta (j,k)\right]}_{T\times q_{\eta }}$, and $\psi ={\left[\psi (k,m)\right]}_{q_{\eta }\times q_{\pi }}$, $j=1,\dots ,T$; $k=1,\dots ,q_{\eta }$; $m=1,\dots ,q_{\pi }$.

Construct the matrix $\widehat{\mathbf{A}}$ in which the entries are $\widehat{A}(k,k^{\prime })=<{\widehat{x}}_{k,p},{\widehat{x}}_{k^{\prime },p}>,$ where ${\widehat{x}}_{k,p}\left(s\right)={\widehat{x}}_k\left({\displaystyle \frac{s+p-1}{p}}\right),0\le s\le 1.$ Define the orthonormal eigenvectors ${\widehat{\beta }}_{\mathbf{k}}$ with corresponding ordered eigenvalues ${\widehat{\xi }}_1\ge \dots \ge {\widehat{\xi }}_{q_{\eta }}$ as $\widehat{\mathbf{A}}{\widehat{\beta }}_{\mathbf{k}}={\widehat{\xi }}_k{\widehat{\beta }}_{\mathbf{k}},1\le k\le q_{\eta },$ and denote $\widehat{\mathbf{B}}=[{\widehat{\beta }}_{\mathbf{1}},\dots ,{\widehat{\beta }}_{{\mathbf{q}}^{*}}]$, where $q^{*}=max\{k\in \{1,\dots ,q_{\eta }\}:\left|\right|{\widehat{z}}_{k,p}{\left|\right|}^2\ge 0.9p\}$ and ${\widehat{z}}_{k,p}\left(s\right)={\sum }_{i=1}^{q_{\eta }}{\widehat{\beta }}_{k,i}{\widehat{x}}_{k,p}\left(s\right)$.

Following [10], the test statistic is constructed as

$${\widehat{\tau }}_p={\displaystyle \frac{1}{T}}{\left(vec\left[{\widehat{\mathbf{B}}}^{\tau }\widehat{\psi }\right]\right)}^{\tau }{\left[({\mathbf{I}}_{q_{\epsilon }}\otimes {\widehat{\mathbf{B}}}^{\tau })(\widehat{\mathbf{C}}\otimes \widehat{\Lambda })({\mathbf{I}}_{q_{\epsilon }}\otimes \widehat{\mathbf{B}})\right]}^{-1}vec\left[{\widehat{\mathbf{B}}}^{\tau }\widehat{\psi }\right],$$

(11)

where $\widehat{\Lambda }=diag({\widehat{\lambda }}_1,\dots ,{\widehat{\lambda }}_{q_{\eta }}),\widehat{\mathbf{C}}={\displaystyle \frac{1}{T}}{(\pi -\eta \widehat{\psi })}^{\tau }(\pi -\eta \widehat{\psi })$. The test statistic ${\widehat{\tau }}_{p+1}$ has an approximately chi-squared distribution with degree $q_{\pi }{q}^{*}$.

4. Theoretical Results

In this section, we discuss the asymptotic properties of both initial and improved estimation of $g_m(u,x_m)$. Moreover, the consistency of estimation of order p is derived. All proofs are given in Appendix.

Throughout the remainder of this article, for any fixed interval $[a,b]$, we denote the space of l-th order smooth functions as $C^{\left(l\right)}[a,b]=\left\{g\right|g^{\left(l\right)}\in [a,b]\}$ and the class of Lipschitz-continuous functions for some fixed constant $C>0$ as $Lip([a,b],C)=\left\{g\right||g\left(x\right)-g\left(x^{\prime }\right)|\le |x-x^{\prime }|,\forall x,x^{\prime }\in [a,b]\}.$ Let $S_{x_m}$ and $S_{z_m}$ denote the supports of $x_m$ and $z_m$, respectively. Then, the supports of $\mathbf{x}$ and $\mathbf{z}$ are $S_x={\prod }_{m=1}^k{S}_{x_m}$ and $S_z={\prod }_{m=1}^k{S}_{z_m}$, respectively. The necessary assumptions for the asymptotic results are as follows.

(A1) For any $d\in \mathcal{F}$, d is differentiable and there exists a constant $M>1$ such that ${\left|\right|d\left|\right|}_{\infty }$, $\left|\right|1/{d\left|\right|}_{\infty }$, and $\left|\right|d^{\prime }{\left|\right|}_{\infty }$ are all bounded by M.

(A2) (a) The kernel density $\mathcal{K}$ is Lipschitz-continuous, bounded, and symmetric about 0. Furthermore, $\mathcal{K}\in Lip([-1,1],L_k)$ for some constant $L_k>0$. (b) The kernel density $\mathcal{K}$ is such that ${\int }_0^1\mathcal{K}\left(u\right)du>0$, ${\int }_{\mathcal{R}}\left|u\right|\mathcal{K}\left(u\right)du$, ${\int }_{\mathcal{R}}{\mathcal{K}}^2\left(u\right)du$, and ${\int }_{\mathcal{R}}\left|u\right|{\mathcal{K}}^2\left(u\right)du$ are finite.

(A3) The covariates $x_{t,m},z_{t,m}$, $1\le m\le k$, and the errors ${\epsilon }_t\left(u\right)$ satisfy the following moment conditions: for some $s>2$,

$$\underset{1\le t\le T}{max}\underset{1\le m\le k}{max}E\left(\right|x_{t,m}{|}^{2s})<\infty ,\underset{1\le t\le T}{max}\underset{1\le m\le k}{max}E\left(\right|z_{t,m}{|}^{2s})<\infty ,\underset{1\le t\le T}{max}\underset{u}{sup}E\left(\right|{\epsilon }_t\left(u\right){|}^{2s})<\infty .$$

For each $t=1,\dots ,T$, the covariance function $Cov({\epsilon }_t\left(s\right),{\epsilon }_t\left(v\right))={\Sigma }_t(s,v)$ has finite nondecreasing eigenvalues ${\lambda }_1\le \dots \le {\lambda }_{max}$ satisfying ${\sum }_j{\lambda }_j<\infty$.

(A4) The additive component functions $g_m(u,x_m),1\le m\le k$, are continuous on $[0,1]\times [a_m,b_m]$ and twice continuously partially differentiable with respect to u and $x_m$, where $[a_m,b_m]$ is a compact subset of $S_{x_m}$.

(A5) $N_0\sim {\left(nT\right)}^{1/6}lognT$, $N_m\sim {\left(nT\right)}^{1/6}lognT$, $1\le m\le k$, $h\sim n^{-1/3}$, as $n,T\to \infty$.

Remark 1.

Assumption (A1) is basic and essential for deriving the consistency of densities after transformation. The conditions in (A2) on the kernel function $\mathcal{K}(·)$ are mild and can be satisfied by commonly used kernel functions such as the uniform and Epanechnikov kernels. The moment conditions in (A3) are crucial for deriving the uniform convergence and other asymptotic properties based on the spline function. The smoothness conditions for the component functions in (A4) are greatly relaxed. The conditions in (A5) are commonly applied in spline smoothing to ensure optimal convergence rates.

We first derive the uniform consistency of the initial estimates of bivariate functions $g_m(u,x_m)$, as stated in Theorem 1.

Theorem 1.

Theorem 1. Assume that assumptions (A1)–(A5) hold, and that ${\tilde{g}}_m(u,x_m)$ are the initial estimates of $g_m(u,x_m)$ defined by (5), $m=1,\dots ,k$. Then, as $n\to \infty$ and $T\to \infty$, it holds that

$$\underset{u,x_m\in [0,1]}{sup}|{\tilde{g}}_m(u,x_m)-g_m(u,x_m)|=O_p\left({\left(nT\right)}^{-1/3}log\left(nT\right)+n^{-1/3}\right).$$

Theorem 2 characterizes the uniform convergence of the improved estimation of $g_m(u,x_m)$, and Theorem 3 describes the asymptotic properties of both the initial and improved estimations.

Theorem 2.

Theorem 2.Assume that assumptions (A1)–(A5) hold, ${\widehat{g}}_m(u,x_m)$ are the improved estimates of $g_m(u,x_m)$ defined by (9), $m=1,\dots ,k$, and that the order of the functional error process p is known. Then, as $n\to \infty$ and $T\to \infty$, the following holds:

$$\underset{u,x_m\in [0,1]}{sup}|{\widehat{g}}_m(u,x_m)-g_m(u,x_m)|=O_p({\left(nT\right)}^{-1/3}{(log\left(nT\right))}^{-2}+n^{-1/3}).$$

To present the asymptotic normality of the estimations, some notations are given as follows. Denote $\mathbf{b}(u,x_{t,m})={(b_{1,1,m}(u,x_{t,m}),\dots ,b_{N_0,N_m,m}(u,x_{t,m}))}^{\tau }$, ${\mathbf{b}}_{\mathbf{z}}(u,x_{t,m})=z_{t,m}\mathbf{b}(u,x_{t,m})$, ${\mathbf{Bz}}_{\mathbf{t},\mathbf{m}}={({\mathbf{b}}_{\mathbf{z}}(u_1,x_{t,m}),\dots ,{\mathbf{b}}_{\mathbf{z}}(u_n,x_{t,m}))}_{n\times N_0{N}_m}^{\tau }$, ${\mathbf{B}}_{\mathbf{m}}={({\mathbf{Bz}}_{\mathbf{1},\mathbf{m}}^{\tau },\dots ,{\mathbf{Bz}}_{\mathbf{T},\mathbf{m}}^{\tau })}^{\tau }$, and $\mathbf{B}=({\mathbf{B}}_{\mathbf{1}},\dots ,{\mathbf{B}}_{\mathbf{k}})$.

Let ${\mathbf{B}}_{*}=\mathbf{B}/\sqrt{nT}$, and ${\mathbf{A}}_{\mathbf{m}}=(\mathbf{0},\dots ,\mathbf{I},\dots ,\mathbf{0})$ is an $1\times k$ block matrix, with the m-th block being $N_0{N}_m\times N_0{N}_m$ identity matrix and the j-th ($j\ne m$) block being $N_0{N}_m\times N_0{N}_j$ zero matrix.

Theorem 3.

Theorem 3.Assume that assumptions (A1)–(A5) hold, ${\tilde{g}}_m(u,x_m)$ and ${\widehat{g}}_m(u,x_m)$ are the initial and improved estimates of $g_m(u,x_m)$ defined by (5) and (9), $m=1,\dots ,k$, respectively. Then, as $n\gg T\to \infty$, the following hold for all $u\in (0,1)$ and $x_m\in [0,1]$:

(i) The initial estimate ${\tilde{g}}_m(u,x_m)$ is asymptotically normally distributed, i.e.,

$$\sqrt{nT}{\left({\mathbf{C}}_{\mathbf{m}}{\Sigma }_{\epsilon }{\mathbf{C}}_{\mathbf{m}}^{\tau }\right)}^{-1}({\tilde{g}}_m(u,x_m)-g_m(u,x_m))\stackrel{D}{\to }N(0,1),$$

where ${\mathbf{C}}_{\mathbf{m}}={\mathbf{b}}^{\tau }(u,x_m)E\left({\mathbf{A}}_{\mathbf{m}}{\left({\mathbf{B}}_{*}^{\tau }{\mathbf{B}}_{*}\right)}^{-\mathbf{1}}{\mathbf{B}}_{*}^{\tau }\right)$, the covariance matrix ${\Sigma }_{\epsilon }={\left({\Sigma }_{\mathbf{t},\mathbf{s}}\right)}_{1\le t,s\le T}$, with ${\Sigma }_{\mathbf{t},\mathbf{s}}=Cov({\epsilon }_{\mathbf{t}},{\epsilon }_{\mathbf{s}})$.

(ii) The improved estimate ${\widehat{g}}_m(u,x_m)$ is asymptotically normally distributed, i.e.,

$$\sqrt{nT}{\left({\mathbf{C}}_{\mathbf{m}}{\Xi }_{\epsilon }{\mathbf{C}}_{\mathbf{m}}^{\tau }\right)}^{-1}({\widehat{g}}_m(u,x_m)-g_m(u,x_m))\stackrel{D}{\to }N(0,1),$$

where the covariance matrix ${\Xi }_{\epsilon }=\mathrm{diag}{\left({\Xi }_{\mathbf{t},\mathbf{t}}\right)}_{1\le t\le T}$, with ${\Xi }_{\mathbf{t},\mathbf{t}}(u,s)={\sigma }_t^2(u,s)$.

5. Numerical Study

In this section, we conduct two simulation studies to demonstrate the performance of the proposed identification and estimation procedure for the additive model.

5.1. Case 1

With the assumption that the auto-regressive error process is known, this case is conducted to demonstrate the performance of the estimation procedure with finite n and T. We consider the following DVCA-FAR(1) model:

$$f_t\left(u\right)=z_{t,1}{g}_1(u,x_{t,1})+z_{t,2}{g}_2(u,x_{t,2})+{\epsilon }_t\left(u\right),0\le u\le 1,$$

(12)

and the error function ${\epsilon }_t\left(u\right)$ is

$${\epsilon }_t\left(u\right)=\int {\gamma }_1(s,u){\epsilon }_{t-1}\left(s\right)ds+e_t\left(u\right),2\le t\le T.$$

Let the bivariate varying-coefficient functions be

$$g_1(u,x_{t,1})=sin\left(2\pi u\right)(2x_{t,1}-1),g_2(u,x_{t,2})=sin\left(2\pi u\right)sin\left(2\pi x_{t,2}\right),$$

and the coefficient functions be

$${\gamma }_1(s,u)=0.2us,e_t\left(u\right)=0.2{\eta }_{t,1}sin\left(\pi u\right)+{\eta }_{t,2}sin\left(2\pi u\right),$$

where ${\eta }_{t,1}\sim N(0,0.1^2)$, ${\eta }_{t,2}\sim N(0,0.{05}^2)$, and ${\eta }_{t,1}$ are independent of ${\eta }_{t,2}$ for $u\in [0,1]$.

The covariates $z_{t,1},z_{t,2}$ are generated by $N(0,1)$ and $N(0,0.5^2)$, respectively, while $x_{t,1},x_{t,2}$ are generated by ${(x_{t,1},x_{t,2})}^{\tau }={(\Phi \left(v_{t,1}\right),\Phi \left(v_{t,2}\right))}^{\tau }$, $1\le t\le T$, where $\Phi$ is the cumulative distribution function of the standard normal distribution and $v_{t,1}$, $v_{t,2}$ are mutually independent variables of the standard normal distribution.

To generate the response densities, for each given $\mathbf{Z}=\mathbf{z}$ and $\mathbf{X}=\mathbf{x}$, let $\alpha (u,\mathbf{x},\mathbf{z})$ be the additive surface given by $\alpha (u,\mathbf{x},\mathbf{z})={\sum }_{m=1}^p{z}_m{g}_m(u,x_m)$. The conditional quantile function $Q(·|\mathbf{x},\mathbf{z})$ with the error process $\epsilon \left(u\right)$ satisfy $Q\left(u\right|\mathbf{x},\mathbf{z})=F^{-1}\left(u\right|\mathbf{x},\mathbf{z})=\theta {(\mathbf{x},\mathbf{z})}^{-1}{\int }_0^{u}exp\{\alpha (v,\mathbf{x},\mathbf{z})+\epsilon \left(v\right)\}dv$, where $\theta (\mathbf{x},\mathbf{z})={\int }_0^{1}exp\{\alpha (v,\mathbf{x},\mathbf{z})+\epsilon \left(v\right)\}dv$.

Implementing the conditional quantile function to $\{U_{t,1}\le \dots \le U_{t,n_t}\}\sim U(0,1)$, which are independent of ${\mathbf{X}}_{\mathbf{t}}$ and ${\mathbf{Z}}_{\mathbf{t}}$, we can get the random samples ${\mathbf{Y}}_{\mathbf{t}}=\{Y_{t,j}=Q\left(U_{t,j}\right|{\mathbf{X}}_{\mathbf{t}},{\mathbf{Z}}_{\mathbf{t}}):1\le j\le n_t\}$ for each $1\le t\le T$, so that $Y_{t,1},\dots ,Y_{t,n_t}\sim d_t$, where $d_t$ is the random response density. Denoting $f_t\left(u\right)=\Psi \left(d_t\left(u\right)\right)$, we obtain the transformed density expressed in model (2). Without loss of generality, we assume that $n_t=n$ independent and identically distributed observations are available for each response distribution.

With $T=100,n=100$, the estimation is conducted with 200 Monte Carlo runs. Figure 2 displays the true curves of the FAR error process $\epsilon \left(u\right)$ in panel (a) and its corresponding spline-based estimations in panel (b). Figure 3 presents a comprehensive view of the true surface of $g_m(u,x_m)$ alongside the average estimations obtained from 200 Monte Carlo simulations. Specifically, the left panel displays the true densities, the middle panel shows the initial spline estimates derived without accounting for the $FAR\left(1\right)$ error process, and the right panel provides the improved estimations achieved through the same method after removing the estimated error process. To better illustrate the performance, the bivariate function estimations are presented from two different relative perspectives, offering a clearer understanding of the improvements made through the error correction process.

Figure 2 indicates that the estimation of $\epsilon \left(u\right)$ yields highly accurate results, which is further corroborated by the findings presented in Figure 3. Specifically, the right panel of Figure 3 illustrates that the improved estimation, after accounting for the FAR error process, significantly outperforms the initial estimation shown in the middle panel. This improvement highlights the effectiveness of the proposed methodology in refining the bivariate function estimates by appropriately addressing the error structure.

For further comparison, we conducted simulations with sample sizes of $T=50,100$ and $n=50,100$ observations. We use the root mean squared errors (RMSEs) to measure the accuracy of the estimations, including the initial and improved estimations of $g_m(u,x_m)$. The RMSE is defined as

$$RMSE\left({\dot{g}}_m\right)={\displaystyle \frac{1}{T}}\sum _{t=1}^T\{{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left|\right|{\dot{g}}_m(u_i,x_{t,m})-g_m(u_i,x_{t,m}){\left|\right|}_2^2{\}}^{{\displaystyle \frac{1}{2}}}.$$

where ${\dot{g}}_m$ is ${\tilde{g}}_m$ or ${\widehat{g}}_m$.

Based on the results from 200 Monte Carlo simulations, Table 1 presents the average root mean square errors (RMSEs) along with their standard deviations for both the initial and improved estimations of $g_m(u,x_m)$. The findings reveal that the RMSEs of the bivariate functions decrease as both the sample size T and the number of observations n increase. Notably, the RMSEs associated with the improved estimation are consistently smaller than those for the initial estimation. This is to be expected, as the initial estimates were derived without accounting for the error process, which, as demonstrated, has a significant impact on the accuracy of the results. By incorporating the error process and removing its estimated effects in the improved estimation, the model yields more accurate and refined results.

To provide a clearer explanation of the theoretical validity differences between the two estimations, we calculate their biases and standard deviations and compare them. Taking the first example in the numerical simulation section as an example.

Table 2 presents the average bias and standard deviation (SD) of both the initial and improved estimations of $g_m(u,x_m)$. The results clearly indicate that, across all settings of the model sample size, the standard deviation of the improved estimates of the bivariate additive functions is substantially smaller than that of the initial estimates, with the bias also being correspondingly reduced. Moreover, as the sample size increases, both the standard deviation and bias of the estimates decrease, further reinforcing the reliability of the improved method. This finding numerically substantiates the claim that the improved estimation method results in a smaller asymptotic variance-covariance matrix compared to the initial estimation, thereby enhancing the precision and robustness of the estimates.

5.2. Case 2

Case 2 is conducted to demonstrate the efficiency of identifying the auto-regressive order of the functional error process. The densities are also generated from model (12), but with the $FAR\left(2\right)$ error function, where ${\gamma }_2(s,u)={\displaystyle \frac{1}{4}}u{s}^2$. All other settings are the same as for Case 1.

Table 3 presents the empirical power of the testing algorithm used to determine the order of the FAR error process under various settings and significance levels. The results clearly demonstrate that the power of the test increases as the sample size T and the number of observations n grow. Specifically, the power approaches 1 as both T and n increase to 100 when testing the null hypothesis of an independent and identically distributed (i.i.d.) sample. This suggests that the test becomes increasingly reliable with larger sample sizes. However, the power is slightly lower when testing the null hypothesis of order 1 against order 2, which is expected given the complexity of distinguishing between these two orders. Furthermore, the size of the test remains low when testing order 2 against a higher-order process, confirming the accuracy and feasibility of the testing algorithm for determining the appropriate order of the functional error process. These findings validate the effectiveness of the proposed testing algorithm in practical applications, ensuring its robustness and precision in various settings.

Furthermore, to assess the efficiency of the auto-regressive order p on the overall estimation results, Table 4 presents the average RMSEs for the bivariate varying-coefficient functions. The observed pattern closely mirrors the results from Case 1, further reinforcing the effectiveness and reliability of the proposed model’s identification and estimation procedures. These findings provide strong empirical evidence that the auto-regressive order plays a crucial role in enhancing the accuracy of the estimation process. The consistency of the RMSE results across different settings underscores the robustness of the model, confirming its ability to effectively account for the error structure and yield precise estimates of the bivariate varying-coefficient functions.

6. Real Data Analysis

In this section, we demonstrate the feasibility and efficiency of the proposed estimation procedure through the analysis of two real-world datasets. By applying the developed methodology to empirical data, we aim to showcase the practical utility of the model in capturing the underlying patterns and dependencies inherent in the datasets. The analysis serves not only to validate the effectiveness of the proposed estimation approach but also to highlight its applicability across different domains, offering insights into the model’s versatility and robustness in handling complex, time-dependent, and non-Euclidean data structures.

6.1. COVID-19 Data

On March 11, 2020, the World Health Organization (WHO) declared COVID-19, an infectious disease caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a global pandemic. The rapid and widespread transmission of the virus led to unprecedented global health challenges, with countries around the world instituting lockdowns and other measures to curb the spread of the disease. As of August 15, 2021, official statistics from the WHO reported a staggering 221,885,822 confirmed cases and 4,583,539 deaths across nearly all countries, reflecting the profound and far-reaching impact of the pandemic. Given the magnitude of this crisis, it is essential for international health organizations and research institutions to closely monitor the evolving global trends of COVID-19. Such monitoring facilitates timely and accurate analysis, enabling effective public health responses and the development of strategies for medical treatment, prevention, and control of future outbreaks. Understanding the dynamics of the epidemic through data-driven models is thus critical for informing policy decisions and improving global health outcomes in the face of such a devastating crisis.

To demonstrate this proposition, we select the mortality rate as an appropriate metric for measuring the global trend of the COVID-19 pandemic. The mortality rate is defined as the ratio of cumulative deaths per day to the total population of each country, serving as a key indicator of the disease’s lethality and spread. Notably, the data required to calculate the mortality rate are inherently temporally dependent, as they rely on the data from the previous day. This results in the mortality rate, and thus the global trend of the epidemic, exhibiting temporal auto-correlation.

The data on COVID-19-related deaths, which are critical for our analysis, can be accessed from the Johns Hopkins University repository. This repository includes a dynamic tracking map that provides a comprehensive view of the global trends related to the pandemic. The dataset, which is publicly available at https://www.jhu.edu/, covers the period from January 22, 2020, to April 15, 2021. Additionally, the most recent total population data for each country, necessary for calculating the mortality rate, can be obtained from the World Bank’s online platform, accessible at https://data.worldbank.org/indicator. These publicly available datasets serve as a rich resource for tracking the progression of the pandemic and for conducting rigorous statistical and epidemiological analyses aimed at understanding the disease dynamics across different regions.

Due to the varying outbreak times across different countries and regions, we define the origin of the time scale as the point at which the cumulative number of confirmed COVID-19 cases reached 100 in all countries. For this analysis, we focus on a dataset containing daily cumulative deaths from 189 countries, considering a 100-day period following this reference time. At each time point t, we estimate the density function of the mortality rate, denoted as ${\widehat{d}}_t\left(y\right)$, using the observations from the 189 countries. Figure 1 (a) presents the estimated densities of the global mortality rate (‰) over the 100-day interval, with data records from up to 189 countries for each time point. Figure 1 (b) offers an alternative perspective by displaying the estimated densities for three selected days. From the figures, it is evident that the mortality rate densities are well-defined across the observed period. Moreover, a temporal dependency among the distributions is clearly observable, which suggests the presence of an auto-regressive process in the data, potentially supporting the hypothesis of a FAR error structure.

The primary goal of this analysis, based on the COVID-19 data, is to identify the FAR process underlying the mortality rate and estimate its component functions. For the sake of simplicity, we begin by considering a special case in which the covariate z is constant (=1) and x represents a time scale, denoted as $t/T$, namely, we consider the following model:

$${\widehat{f}}_t\left(u\right)=\Psi \left({\widehat{d}}_t\right)\left(u\right)=g_1(u,x_{t,1})+{\epsilon }_t\left(u\right),1\le t\le 100,$$

where ${\epsilon }_t\left(u\right)={\sum }_{l=1}^p\int {\gamma }_l(u,s){\epsilon }_{t-l}\left(s\right)ds+e_t\left(u\right)$ and $x_{t,1}$ denotes the time scale $t/T$.

Based on the initial spline estimation of $g_1(u,x_1)$, the testing algorithm is employed to determine the order of the FAR process. Table 5 presents the test results, specifically the p-values under different hypotheses. The table reveals that the observed p-values indicate significant evidence of auto-correlation in the data. This suggests that the underlying process can indeed be effectively modeled as a first-order functional auto-regressive process, denoted as $FAR\left(1\right)$. Such findings provide strong empirical support for the presence of temporal dependencies in the COVID-19 mortality rate, further justifying the application of a FAR error structure in modeling the dynamics of the epidemic.

Figure 4 presents the heat map of the estimated bivariate function $g_1(u,x_1)$ after accounting for the functional error process and identifying the auto-regressive order. The heat map reveals a relatively stable temporal pattern, with the function initially reaching a minimum value at lower values of u, gradually increasing, and eventually reaching a maximum at later time points. This pattern highlights the underlying dynamics of the COVID-19 mortality rate across successive days. The observed correlation between the data from consecutive days supports the notion that the global mortality rate exhibits significant temporal dependencies. This is consistent with the design of the mortality rate measure, which is derived from previous daily data, reflecting the evolving trend of the pandemic on a global scale. The results further validate the presence of auto-correlation in the mortality rate data and underscore the necessity of incorporating a functional auto-regressive process to capture these dependencies effectively.

6.2. USA Income Data

Personal income statistics play a crucial role in enabling governments to understand the dynamics between national income, spending, and saving, while also serving as an important tool for evaluating and comparing the economic well-being across different regions or nations. In this context, we focus on the density time series of per capita personal income, which is defined as the total personal income of an area divided by its population. This measure provides a more granular perspective on the economic conditions within a region, reflecting the distribution and trends in income on a per-person basis over time. By analyzing such time series, policymakers and researchers can gain insights into the long-term economic trajectory of a region, assess disparities in income distribution, and make informed decisions regarding fiscal policies, social welfare programs, and economic development strategies.

Income data for the USA are publicly available at the official website of the United States Bureau of Economic Analysis (http://www.bea.gov/). We consider the quarterly per capita personal income of 50 states in the USA from the first quarter of 2010 to the fourth quarter of 2020, namely, $t=1,\dots ,44$. For each t, we obtain the density function of per capita income, ${\widehat{d}}_t\left(y\right)$, each with 50 observations. As the quarterly personal income is an economic measure based on national conditions, we choose two related covariates, `GDP’ (quarterly gross domestic product of the USA) and `Population’ (quarterly total population of the USA), which can also be obtained from http://www.bea.gov/.

The income curve, traditionally studied in economics as panel data, primarily reflects the relationship between consumers’ equilibrium points. As individuals’ income levels fluctuate, the connections between these equilibrium points form a trajectory, symbolizing not only an increase in income but also a corresponding rise in consumer satisfaction. This approach emphasizes the dynamic nature of income growth and its impact on individual well-being, providing valuable insights into consumer behavior over time.

In contrast, the income density curve, treated as functional data, is used to describe the distribution of income within a specific region or demographic group. It offers a graphical representation of the characteristics and trends in income across various income intervals, providing a more holistic view of the socio-economic landscape. By analyzing the income density curve, the degree of income inequality within a population can be effectively observed, revealing important patterns of wealth distribution. This type of curve is crucial for economic research as it facilitates a deeper understanding of consumption behavior, socio-economic conditions, and the formulation of societal policies.

Moreover, income density curves are instrumental in economic forecasting and analysis. By examining trends in income distribution, economists can combine insights into consumer preferences and consumption habits at different income levels. This integration enables predictions about future economic conditions and shifts in consumption patterns, making the income density curve a key tool for anticipating changes in both microeconomic and macroeconomic contexts. In this way, the income density curve plays an essential role in informing policy decisions, economic strategies, and the broader understanding of economic well-being.‌

Figure 5 (a) illustrates the density time series of quarterly personal income over a period of 44 quarters. The density curves reveal that, over the past decade, the overall distribution of per capita income across various states in the United States has exhibited a consistent pattern. Specifically, there are relatively few individuals in the high-income and middle-to-high-income brackets, a moderate number in the middle-income category, and a larger proportion in the middle-to-low-income segments.

To further elucidate the temporal evolution of income distribution, Figure 5 (b) presents the density curves for three distinct time points: the second quarter of 2015, the first quarter of 2017, and the third quarter of 2018. A clear trend emerges, showing a gradual shift towards higher income values over time, accompanied by a corresponding decrease in the peak of the density curve. This shift is unsurprising given the broader economic and technological advancements in modern society. As the economy progresses, the proportion of low-income individuals in the United States has been steadily declining, while the number of middle-to-high-income individuals has been rising. Consequently, the distribution of income is becoming more balanced, with a growing proportion of the population occupying the middle to high-income brackets. This pattern reflects the broader trends of economic development and income redistribution that have taken place in recent years.

We consider the following DVCA-FAR model:

$${\widehat{f}}_t\left(u\right)=\Psi \left({\widehat{d}}_t\right)\left(u\right)=g_0(u,x_{t,0})+z_{t,1}{g}_1(u,x_{t,1})+z_{t,2}{g}_2(u,x_{t,2})+{\epsilon }_t\left(u\right),1\le t\le 44,$$

where ${\epsilon }_t\left(u\right)={\sum }_{l=1}^p\int {\gamma }_l(u,s){\epsilon }_{t-l}\left(s\right)ds+e_t\left(u\right)$. $z_{t,1}$ denotes the quarterly gross domestic product of the USA, $z_{t,2}$ denotes the quarterly total population of the USA, and $x_{t,0},x_{t,1},x_{t,2}$ denote the time scale $t/T$.

Similar to the previous estimation procedure, the testing algorithm is applied to determine the optimal order based on the initial spline estimations. Table 6 presents the p-values of the test under various hypotheses, which indicate that the $FAR\left(2\right)$ model is the most appropriate for modeling the error process in this context. This result suggests that a second-order functional autoregressive process effectively captures the autocorrelated structure of the error terms in the income data.

Using the three-step estimation procedure, estimations of the bivariate varying-coefficient functions can be obtained. Figure 6 displays the heat maps of the three bivariate functions, where $g_0(u,x_0)$ represents the common effect of the data over time, $g_1(u,x_1)$ and $g_2(u,x_2)$ reflect the impact of quarterly GDP and quarterly total population on data respectively. The heat map of $g_0$ reveals a clear pattern, where high and low values alternate over time, indicating that individuals with both higher and lower per capita income experience similar effects, in contrast to those in the middle-income bracket, who show an opposing trend. For $g_1$, the figure demonstrates a mode that is consistent for both small and large values of u, but varies over time. Initially, the function reaches a maximum, then decreases to a minimum before increasing again towards the end of the time scale. In contrast, the effect of population, as shown in the heat map of $g_2$, exhibits a trend that is opposite to that of the common effect. The population impact is relatively consistent across both higher and lower per capita income groups. Taken together, these results suggest a significant dependency of the quarterly personal income distributions in the United States on previous statistics, with the dynamics of income closely tied to both macroeconomic factors (such as GDP) and demographic factors (such as population).

7. Discussion

Data collected from sequential time points often exhibit autocorrelation, which must be addressed in the modeling process. Additionally, the analysis of non-Euclidean data has become increasingly common. To tackle these challenges, we propose a varying-coefficient additive model with density responses, incorporating a Functional Auto-Regressive (FAR) error process. Due to the complexity and constraints of the data, we first map density functions into a linear space using the transformation method by [12]. We then develop a three-step estimation procedure for the varying-coefficient components. First, B-spline series approximation is used to obtain initial estimates of the bivariate varying-coefficient components without considering the functional errors. Next, the FAR error process order is determined using the test statistic from [10] based on the residuals from the initial estimation. Finally, the FAR error process is removed, and improved spline estimators are constructed for the varying-coefficient components. Theoretical results, including convergence rates and asymptotic properties, are derived for both the initial and improved estimations. The performance of the proposed method is demonstrated through simulations and two real datasets, showing the importance of accounting for autocorrelation and validating the efficiency of the proposed approach.

Future research can explore a variety of related problems. In this study, we employed the varying-coefficient additive model to establish the relationship between density function responses and scalar predictors. With the advent of large and complex datasets, functional predictors have become increasingly valuable in analyzing practical data applications. Furthermore, to account for sequence dependence, we utilized the FAR process to model the correlation structures in the data. However, more intricate models representing autocorrelation could also be explored to further refine the representation of sequence dependence in functional errors. Future work will focus on extending these approaches and conducting additional studies in these areas.