Forecasting the monthly retail sales of electricity based on the semi-functional linear model with autoregressive errors

1.   Introduction

Functional regression modeling has been extensively studied as an important research field in functional data analysis (FDA). For example, to investigate the linear relationship between a scalar response and a functional explanatory variable, comprehensive studies focused on functional linear regression models have been conducted by Manteiga and Calvo ^[1] and Bande et al. ^[2]. However, the linear assumption is sometimes too restrictive to describe the data structure comprehensively. To this end, the functional nonparametric regression model, as an extension of the nonparametric model in the functional data setting, has also been proposed (see Ferraty and Vieu ^[3] and Belarbi et al. ^[4]). Furthermore, to improve the power of prediction and interpretation of functional regression models, various functional semiparametric models have been proposed by introducing some additional real-valued predictors. For example, Pérez and Vieu ^[5] proposed a semi-functional partial linear regression model by combining a traditional linear model with a nonparametric treatment of functional data. Shin ^[6] and Shin and Lee ^[7] studied the problem of estimation and prediction of the partial functional linear regression model, respectively. Yao et al. ^[8] extended the partial functional linear regression to quantile regression settings and established a more flexible and robust method. Kong et al. ^[9] investigated a class of partial function linear regression models for processing multiple functional and non-functional explanatory variables and automatically identifying essential risk factors through appropriate regularization. In addition, Zhou and Chen ^[10] first introduced the semi-functional linear model (SFL) by combining the features of a functional linear regression model and a traditional nonparametric model and studied its spline estimation. Peng et al. ^[11] proposed a varying coefficient partially functional linear regression model, which combines a varying coefficient regression model and a functional linear regression model, and investigated the estimation problem based on the polynomial spline method. Novo et al. ^[12] developed an automatic and location-adaptive estimation procedure for the functional single-index model. For a recent survey on this field, we refer to some literature like Cuevas ^[13], Goia and Vieu ^[14], Ling and Vieu ^[15,16], and Aneiros et al. ^[17].

Meanwhile, given that electricity cannot be stored, modeling and forecasting electricity demand and price are of primary interest to electric market agents. Traditional prediction models, like simple time series models, often fail to address the complexities of electricity demand, so functional regression models have been applied to electricity demand forecasting. For example, Vilar et al. ^[18] utilized functional nonparametric and semi-functional partial linear models to forecast next-day electricity demand and price. Based on the robust functional principal component analysis and nonparametric models with both functional response and covariate, Aneiros et al. ^[19] proposed two functional regression models to forecast daily curves of electricity demand and price. In addition, using residual-based bootstrap algorithms, Vilar et al. ^[20] provided two procedures to obtain prediction intervals for electricity demand and price based on the functional nonparametric autoregressive model and the partial linear semiparametric model. Peláez et al. ^[21] also proposed three model-based procedures to construct bootstrap prediction regions for daily curves of electricity demand and price.

All the models mentioned above are based on a critical assumption that the random errors are independently and identically distributed. However, the independence assumption may be inappropriate in many applications. For example, the electricity consumption dataset that has been extensively studied by Pérez and Vieu ^[22], Yu et al. ^[23], and Yang et al. ^[24] exhibits a specific autocorrelation structure. If this autocorrelation structure is ignored, it may lead to inaccurate prediction. To this end, many works have been developed to deal with the autocorrelation structure in functional regression models. For example, Niang and Guillas ^[25] studied parametric estimation of a functional semiparametric partially linear model with autoregressive errors. Zhang et al. ^[26] proposed a functional polynomial regression with autoregressive errors. Yu et al. ^[23] considered composite quantile estimation for the partial functional linear regression model with random errors from a short-range dependent and strictly stationary linear process. Xiao and Wang ^[27] also studied the estimation problem of the partial functional linear model with autoregressive errors. Wang et al. ^[28] studied the multiple functional linear model with autoregressive errors. Although these functional regression models with autoregressive errors have improved prediction accuracy, the linear assumption on the covariates is sometimes inappropriate. Therefore, it is necessary to develop more general models to fit the more complex data structures. To achieve this goal, by relaxing the linear assumption on the partial functional linear model, a semi-functional linear model with autoregressive errors is developed in this paper. Based on the functional principal component analysis and the spline method, we obtain the estimators of the slope function, nonparametric function, and autoregressive coefficients, respectively.

The rest of this paper is organized as follows. Section 2 describes the semi-functional linear model with autoregressive errors and presents the estimation method and algorithm. In Section 3, we study the asymptotic properties of the proposed estimators under some regular conditions. A Monte Carlo simulation is presented in Section 4 to illustrate the finite sample performance. The proposed method is illustrated by electricity consumption data in Section 5. Sections 6 and 7 present future prospects and summarize the paper, respectively. All technical details and proofs are provided in the appendix.

2.   Model and estimation

2.1.   Model and estimation

Suppose that $(X_t, Y_t, Z_t), 1\leq t\leq n$, is a random sample from the semi-functional linear model with autoregressive errors (SFLAR)

where $Y_t$ is a real-valued random response variable and $Z_t$ is a real-valued covariate defined on a compact interval $Z = [\underline{z}, \overline{z}]$. The functional covariate $X_t(s)$ defined on an interval $\mathcal{T} \subset R$ is a zero mean, second-order (i.e., $\mathbb{E}|X(s)|^{2} < \infty$ for all $s \in \mathcal{T}$) stochastic process defined on $({\Omega}, B, P)$ with sample paths in the Hilbert space $L^{2}(\mathcal{T})$. The inner product of $L^{2}(\mathcal{T})$ is defined by $\langle u, v \rangle = \int_\mathcal{T}{u(s)v(s)}ds$, for any $u, v \in L^{2}(\mathcal{T})$ and norm $\|u\| = {\langle u, u \rangle}^{1/2}$. The unknown slope function $\beta(s)$ belongs to $L^{2}(\mathcal{T})$. For simplicity, throughout this paper, we assume $\mathcal{T} = [0, 1]$. The nonparametric function $g(\cdot)$ is an unknown smooth function. $\varepsilon_t$ is a random error with zero mean and finite variance independent of $(X_t, Z_t)$. $e_t$'s are independent and identically distributed with mean 0 and variance $\sigma^2$. Note that model (2.1) includes the partial functional linear model with autoregressive errors studied by Xiao and Wang ^[27].

Since the unknown functions $\beta (\cdot)$ and $g(\cdot)$ are infinite-dimensional parameters, it is impossible to obtain their estimators based on finite observation data. Therefore, we need to apply dimension reduction methods. Specifically, if we denote the covariance function of the functional variable $X$ by $C_X (s, u) = Cov[X(s), X (u)]$, and assume that $C_X (s, u)$ is continuous on the interval $[0, 1]^2$, according to the Mercer theorem, we have

where $\lambda_ 1 > \lambda_2 > \cdots > 0$ is the ordered eigenvalue sequence and $\phi_ 1, \phi_2, \cdots$ is the corresponding orthonormal eigenfunction sequence of the covariance operator with kernel $C_X (s, u)$. By the Karhunen-Loève representation, functional variables $X_t(s)$ and the slope function $\beta (s)$ can be respectively expanded as:

where $\xi_{ti} = \int_{0}^{1}{{X_t}(s)\phi_{i}(s)}ds$ are called the principle component scores satisfying $\mathbb{E}\xi_{ti} = 0$ and $\mathbb{E}\xi_{ti}\xi_{tk} = \lambda_iI(i = k)$, and coefficients $b_{i} = \int_{0}^{1}{\beta(s)\phi_{i}(s)}ds$. $\left\lbrace \phi_{i}(s)\right\rbrace _{i = 1}^{\infty}$ forms a standard orthogonal basis of $L^{2}(\mathcal{T})$, commonly referred to as the functional principal component (FPC) basis. Truncate Eqs (2.2) at $M$, and the model (2.1) can be approximated as

Furthermore, we can approximate the nonparametric function $g(\cdot)$ using polynomial splines. Specifically, let $G_{p, J_n}$ be the space of polynomial splines defined on the interval $[\underline{z}, \overline{z}]$ with degree $p-1$ and knot sequence $\underline{z} < z_{1} < \cdots < z_{J_n} < \overline{z}$. The space $G_{p, J_n}$ is a $J$-dimensional linear space, where $J = J(n) = J_n+p$. Following the arguments of de Boor ^[29], we can conclude that, if it is sufficiently smooth, the nonparametric function $g(Z_t)$ can be approximately expressed as

where $B_j, j = 1, 2, \cdots, J$, are the B-spline basis functions in $G_{p, J_n}$. Substituting the approximation of $g(Z_t)$ into the model (2.3) yields

However, since the covariance function $C_X(s, u)$ is unknown, the eigenfunctions $\phi_i$ and variables $\xi_{ti}$ are also unknown and unobservable. To address this issue, we define the empirical version of the covariance function $C_X(s, u)$ by

By Mercer's theorem, we also have

where $(\hat{\lambda}_{i}, \hat{\phi}_{i})$ are the eigenvalue-eigenfunction pairs for the linear operator with kernel $\widehat{C}_X(s, u)$, ordered such that $\hat{\lambda}_1 \geq \hat{\lambda}_2 \geq \cdots \geq 0$. We consider $(\hat{\lambda}_{i}, \hat{\phi}_{i})$ and $\widehat{\xi}_{ti} = \int_{0}^{1} X_t(s) \widehat{\phi_{i}}(s) ds$ as estimators of $(\lambda_{i}, \phi_{i})$ and $\xi_{ti}$, respectively. Then, the estimators $\hat{\beta}(s) = \sum_{i = 1}^{M}\hat{b}_i\hat{\phi}_{i}(s)$ and $\hat{g}(Z_{t}) = \sum_{j = 1}^{J}\hat{c}_{j} B_{j}(Z_t)$ can be obtained by minimizing the following loss function:

where $\boldsymbol{\widehat{\xi_t}} = (\widehat{\xi}_{t1}, \widehat{\xi}_{t2}, \cdots, \widehat{\xi}_{tM})^{\top}$, $\boldsymbol{B_t} = \left(B_{1}(Z_t), B_{2}(Z_t), \cdots, B_{J}(Z_t)\right)^{\top}$, $\boldsymbol{a} = (a_1, a_2, \cdots, a_q)^{\top}$, $\boldsymbol{b} = (b_1, b_2, \cdots, b_{M})^{\top}$, $\boldsymbol{c} = (c_1, c_2, \cdots, c_{J})^{\top}$, and $t = q+1, q+2, \cdots, n$.

To implement the proposed estimation method, we must choose the degrees of the spline functions, the positions, and the number of knots and eigenfunctions. This paper uses a B-spline basis with equally spaced knots and the fixed degree 2. Then, we only need to choose the numbers of B-spline functions ($J$) and eigenfunctions ($M$). Various methods, such as the Akaike information criterion (AIC) and Bayesian information criterion (BIC), can be used to select the truncation parameters $J$ and $M$. This paper uses the BIC criterion to choose the number of $J$ and $M$. The selection of $J$ and $M$ is determined by minimizing the following BIC criteria:

where $\hat \varepsilon_t = Y_t-\sum_{i = 1}^{M}\hat{b}_{i}\hat{\xi}_{ti}-\sum_{j = 1}^{J}\hat{c}_{j} B_{j}(Z_t).$

2.2.   Algorithm

This subsection mainly introduces the algorithm to estimate the unknown parameters. For ease of description, we rewrite the loss function (2.5) as

where $\boldsymbol{\widehat{U}_t} = (\boldsymbol{\widehat{\xi}_t^{\top}}, \boldsymbol{B_t^{\top})^{\top}}$, $\boldsymbol{\widehat{\xi}_t} = (\widehat{\xi}_{t1}, \widehat{\xi}_{t2}, \cdots, \widehat{\xi}_{tM})^{\top}$. Since there is no closed-form expression of the minimizer of the loss function (2.7), we adopt the two-step iteratively least-squares (TSILS) algorithm to estimate the unknown parameters. Specifically, we first set $\boldsymbol{{a}^{(0)}} = (0, 0, \cdots, 0)^{\top}$ and compute $\boldsymbol{{\theta}^{(0)}}$ by minimizing $L_n(\boldsymbol{a^{(0)}}, \boldsymbol{\theta})$. Afterward, we substitute $\boldsymbol{{\theta}^{(0)}}$ into the loss function (2.7) to get $\boldsymbol{{a}^{(1)}}$ by $L_n(\boldsymbol{a}, \boldsymbol{\theta^{(0)}})$, and substitute $\boldsymbol{{a}^{(1)}}$ into (2.7) to get $\boldsymbol{{\theta}^{(1)}} = \arg \min L_n(\boldsymbol{a^{(1)}}, \boldsymbol{\theta})$. Repeat this process until convergence. The TSILS algorithm is summarized as follows:

Algorithm 1. The TSILS algorithm initialization: Given the initial value: $\boldsymbol{{a}^{(0)}} = (0, 0, \cdots, 0)^{\top}$, compute $\boldsymbol{{\theta}^{(0)}} = ({\boldsymbol{\widehat{U}^{\top}}}{\boldsymbol{\widehat{U}}})^{-1}{\boldsymbol{\widehat{U}^{\top}}}Y$ based on the least squares method, where $Y = (Y_1, Y_2, \cdots, Y_n)^{\top}$ and ${\boldsymbol{\widehat{U}^{\top}}} = ({\boldsymbol{\widehat{U}_1}}, {\boldsymbol{\widehat{U}_2, \cdots, {\boldsymbol{\widehat{U}_n}}}})$.

Step 1. Update $\tilde{\boldsymbol{V}}^{(k)}$, $\tilde{\boldsymbol{H}}^{(k)}$:

where $\hat{\varepsilon}_t^{(k)} = Y_t-\boldsymbol{\widehat{U}^{\top} _t}\boldsymbol{\hat{\theta}^{(k-1)}}$, $t = 1, 2, \cdots, n.$

Step 2. Compute ${\boldsymbol{a}}^{(k)}$:

Step 3. Update ${V}^{(k)}_t$, $\boldsymbol{H}_t^{(k)}$:

Step 4. Compute ${\boldsymbol\theta}^{(k)}$:

where

Step 5. Set $k = k+1$, and repeat the above steps 1–4 until convergence. Then, the estimation of the parameter vectors $\boldsymbol{\hat{a}} = (\hat{a}_1, \hat{a}_2, \cdots, \hat{a}_q)^{\top}$, $\boldsymbol{\hat{b}} = (\hat{b}_1, \hat{b}_2, \cdots, \hat{b}_{M})^{\top}$, and $\boldsymbol{\hat{c}} = (\hat{c}_1, \hat{c}_2, \cdots, \hat{c}_{J})^{\top}$ can be obtained.

3.   Theoretical properties

In this section, we investigate the theoretical properties of the proposed estimators. First, we assume that the true slope function, nonparametric function, and autoregressive coefficient of the model are $\beta^{0}(t), g^{0}(z)$, and $a^{0} = \left(a_{1}^{0}, a_{2}^{0}, \cdots, a_{q}^{0}\right)^{\top}$, respectively. For convenience, the following notation is needed. For two sequences of positive numbers $a_n$ and $b_n$, $a_n\lesssim b_n$ signifies that $a_n/b_n$ is uniformly bounded and $a_n\asymp b_n$ if $a_n\lesssim b_n$ and $b_n\lesssim a_n$. Throughout this paper, the constant $C_i$ may change from line to line. In order to establish the theoretical properties of the estimators, the following assumptions need to be introduced:

Assumption 1. All roots of the polynomial $A(z) = 1-\sum_{l = 1}^{q} a_{l} z^{l}$ lie outside the unit circle in the complex plane.

Assumption 2. $\left\{X_{t}, Z_{t}\right\}_{t = 1}^{n}$ and $\left\{e_{t}\right\}_{t = 1}^{n}$ are independent random sequences and $\left\{X_{t}, e_{t}\right\}_{t = 1}^{n}$ has the finite fourth-order moment, respectively, i.e. $\mathbb{E} \left(\int X_{t}^{4}(s)ds\right) < \infty$, $\mathbb{E} \left(e_{t}^{4}\right) < \infty$, for any $1 \leq t \leq n$, and random variable $\xi_{ti}$ satisfy $\mathbb{E}\left(\xi_{ti}^{4}\right) \leq C_1 \lambda_{i}^{2}, i\geq 1$. The density function $f_Z(z)$ of the random variable $Z$ is a continuous bounded function defined on the interval $[\underline{z}, \overline{z}]$. Meanwhile, the random variable $e_t$ is independent of $\left(X_t, Z_t\right)$ and $\mathbb{E} \left(X|Z \right) = 0$.

Assumption 3. The eigenvalues ${\lambda_i}(i = 1, 2, \cdots)$ of the covariance function $C_X$ and the score coefficients $b_i$ satisfy the following conditions, respectively.

$\mathbf{(1)}$ There exist some strictly positive constants $C_2$ and $\alpha > 1$ such that $C_2^{-1} i^{-\alpha} \leq \lambda_{i} \leq C_2 i^{-\alpha}$ and $\lambda_{i}-\lambda_{i+1} \geq C_2 i^{-\alpha-1}, i\geq 1$.

$\mathbf{(2)}$ There exist some strictly positive constants $C_3$ and $\gamma > \alpha/2 + 3$ such that $\left|b_{i}\right| \leq C_3 i^{-\gamma}, i\geq 1.$

Assumption 4. Let $0 < \nu \leq 1$. The nonparametric function $g(z)$ is a continuous differentiable function of order $r$, and for $z_1, z_2 \in [\underline{z}, \overline{z}]$, there are $\left|g^{(r)}(z_1)-g^{(r)}(z_2)\right| \leq C_4|z_1-z_2|^{\nu}$. Define $p = r+\nu$ as a smoothness measure of the function $g(z)$, which satisfies $p \geq \gamma+ (\alpha-1)/2$.

Assumption 5. Define the matrix

where $U_t = (\xi_t^{\top}, B_t^{\top})^{\top}$, $\xi_t = (\xi_{t1}, \xi_{t2}, \cdots, \xi_{tM})^{\top}$, $B_t = (B_{1}(Z_t), B_{2}(Z_t), \cdots, B_{J}(Z_t))^{\top}$, and the matrix $V$ is a positive definite matrix. Denote the minimum eigenvalue of $\Lambda$ by $\rho_{\min }(\Lambda)$, which satisfies $\rho_{\min}(\Lambda) = O_{p}\left(M^{-\alpha}\right)$.

Assumption 6. The tuning parameter $M$, $J$ satisfies $M \asymp J \asymp n^{1/(\alpha+2\gamma)}$.

Remark 1. Assumption 1 is an invertibility condition on the autoregression process and ensures ${\varepsilon_t}$ is a stationary series. Assumption 2 is standard in the literature about the least squares method, which ensures the consistency of $C_ X (t, s)$. Assumption 3 is established for the consistency of the estimated principal component scores. Assumption 3(1) prevents the spacings among eigenvalues from being too small for identification of the slope function $b_i$. Assumption 3(2) makes the slope function sufficiently smooth relative to the covariance function $C_ X (t, s)$. Assumption 4 ensures that the nonparametric function $g (z)$ is sufficiently smooth so that the spline function can approach it. Assumption 5 is a necessary condition for the model consistency to be established. Assumption 6 gives the order of the truncation parameters $M$ and $J$ to obtain the convergence speed of $\hat{\beta}(s)$ and $\hat{g}(z)$, and Assumption 6 requires that the tuning parameters $M$ and $J$ do not grow too fast.

Theorem 1. Under the Assumptions 1–6, we have:

Remark 2. The theorem shows that the convergence rate of the autoregressive coefficients cannot reach the optimal convergence rate $O_{p}\left(n^{-1/2}\right)$ (such as Wang et al. ^[30]), which may be due to the existence of functional infinite dimensional explanatory variables. The convergence rates of $\hat{\beta}(s)$ and $\hat{g}(z)$ are the same as that of the case where the random error is independent and identically distributed (such as Tang ^[31]), indicating that the autocorrelation of the random error does not affect the convergence rate of the slope function and the nonparametric function.

4.   Simulation study

In this section, we conduct a simulation under different settings to investigate the finite sample performance of our proposed method. The simulation data $\{(X_t, Z_t, Y_t), 1\leq t \leq n\}$ are generated from the following SFLAR model:

where $Z_t$ is distributed uniformly on [0, 1], and the functional predictor $X_t = \sum_{j = 1}^{50}{\xi_{tj}\phi_{j}}$, where $\xi_{tj}$ are distributed as independent normal with a mean 0 and variance $\lambda_j = ((j-0.5)\pi)^{-2}$ and $\phi_j = \sqrt{2}sin((j-0.5)\pi t)$. For the unknown functions $\beta(\cdot)$ and $g(\cdot)$, and autoregressive coefficients $(a_1, \cdots, a_q)$, two scenarios have been considered as follows:

Scenario 1. $\beta(u) = \sqrt{2}sin(0.5\pi u)+3\sqrt{2}sin(1.5 \pi u)$, $g(z) = sin(2\pi z)$, $\varepsilon_t = a \varepsilon_{t-1}+e_t$, and the autoregressive coefficient $a$ takes each of the values in the set (–0.5, 0.1, 0.8).

Scenario 2. $\beta(u) = \sqrt{2}sin(0.5 \pi u)+2\sqrt{2}sin(1.5\pi u)$, $g(z) = cos(2\pi z)$, $\varepsilon_t = a_1\varepsilon_{t-1}+a_2 \varepsilon_{t-2}+e_t$, and the autoregressive coefficients $(a_1, a_2)$ take each of the values in the set $\{ (0.1, 0.8), (0.5, -0.8)\}$.

The error variables $e_t$ are $N(0, \sigma^2)$. In both scenarios, $\sigma$ takes 0.5 and 1, corresponding to different signal-to-noise ratios. To implement the proposed method, we apply the BIC criterion of Eq (2.6) to choose the number $J$ of B-spline functions and the number $M$ of eigenfunctions. Specifically, $J$ and $M$ are obtained by selecting the minimum BIC value within the preset range. In this simulation, $J$ is chosen between 3 and 8, and $M$ is selected between 1 and 6. In addition, the errors of the estimators of function $\beta(u)$, $g(z)$ and the model prediction error were evaluated by the integral of square errors $\text{ISE}_1 = \int_0^1 [\hat \beta(u)-\beta (u)]^2 du$, $\text{ISE}_2 = \int_0^1 [\hat g(z)-g (z)]^2 dz$ and the mean square prediction error (MSPE): $\mathrm{MSPE} = \frac{1}{n} \sum_{t = q+1}^{n}\left(\widehat{Y_{t}}-g(Z_{t})-\int_0^1 \beta(u) X_{t}(u) \mathrm{d} u\right)^{2}$, respectively.

In our simulation, we considered sample sizes of $n = 400$, $800$, and $1200$, with each setting repeated 500 times. Tables 1 and 2 summarize the performance of various estimators regarding bias (Bias), root mean square error (RMSE) for the estimated $a$ and $\sigma$, and integrated square errors (ISE) for the estimated $\beta(\cdot)$ and $g(\cdot)$ under normal error conditions. From these tables, we can draw the following observations:

(1) Given the sample size, as the value of $\sigma$ increases, the bias (Bias) and root mean square error (RMSE) of $\hat{a}$ and $\hat{\sigma}$ increase, and the mean integral square error (ISE) of the functions $\hat\beta(u)$ and $\hat g(z)$ increase, indicating that with the rise of variance, that is, the signal-to-noise ratio decreases, the estimation accuracy of the model parameters will decrease. The mean square prediction error (MSPE) of the model will increase.

(2) As the sample size $n$ increases from $n = 400$ to $800$ and $1200$, we see a decrease in the root mean square error (RMSE) and the meawe see n integral square error (ISE), and the performance of the estimation method becomes better and better. When the sample size $n$ increases from $n = 400$ to $800$ and $1200$, the estimation of the autoregressive coefficient is asymptotically unbiased. The results show that the estimation of the autoregressive coefficient $a$, and the model parameters $\hat\beta(u)$ and $\hat g(z)$, are consistent estimators. The estimation curves for the slope and nonparametric function of the model in Scenario 1 ($AR(1)$) with parameters $n = 400$, $a = -0.5$, and $\sigma = 0.5$, and in Scenario 2 ($AR(2)$) with $n = 400$, $a = (0.5, -0.8)$, and $\sigma = 0.5$, were generated from 500 repetitions (see Figures 1 and 2). In the top left of Figures 1 and 2, the true slope function curve $\beta(\cdot)$ (red) and the median estimation curve $\hat{\beta}(\cdot)$ (blue) are displayed for Scenarios 1 and 2, respectively. The top right figure shows a comparison between the true curve (red) and the median estimation curve $\hat{g}(\cdot)$ (blue) of the nonparametric function $g(\cdot)$ in both scenarios. The bottom of Figures 1 and 2 displays the comparison of the estimation curves for the slope and the nonparametric function, repeated 500 times in Scenarios 1 and 2, respectively. Figures 1 and 2 show that when $n = 400$ and $\sigma = 0.5$, the slope $\beta(\cdot)$, the nonparametric function $g(\cdot)$, and the median estimation curve (in blue) closely align with the true curve (in red) for both $AR(1)$ and $AR(2)$. For $\sigma = 0.5$, $a = -0.5$, and $a = (0.5, -0.8)$, the 95% pointwise confidence bands for $\beta(\cdot)$ and $g(\cdot)$ are displayed for sample sizes $n = 400$, $800$, and $1200$ in Scenarios 1 and 2, respectively.

When $n = 400$, $800$ and $1200$, and $\sigma = 0.5$, Figures 3 and 4 display the 95% pointwise confidence bands for Scenarios 1 and 2, respectively. Regardless of the changes in the model parameters $\beta(u)$ and $g(z)$ or the errors from $AR(1)$ or $AR(2)$, the proposed estimation method yields estimates for $\beta(u)$ and $g(z)$ that closely align with the true curve. Moreover, as the sample size increases, $\hat{\beta}(u)$ and $\hat{g}(z)$ remain consistent estimates, demonstrating that the proposed model and estimation method perform well in finite samples.

According to the suggestions of an anonymous reviewer, we run a new simulation to compare the predictive ability of the proposed model with that of existing models. Specifically, the simulation data $\{(X_t, Z_t, Y_t), 1\leq t \leq n\}$ are generated according to the following two scenarios:

Scenario 3. It is the same as Scenario 1 with autoregressive coefficient $a = 0.5$ and standard deviation $\sigma = 0.5$.

Scenario 4. It is the same as Scenario 3 except for $g(Z) = 2Z$.

For each simulation case, we set the sample size to $n = 400$, and the simulation data is divided into the training data $(X_t, Z_t, Y_t)_{t = 1}^{320}$ and the test data $(X_t, Z_t, Y_t)_{t = 321}^{400}$. Similar to Aneiros and Vieu ^[22] and Yu et al. ^[23], we utilize the following two error criteria: mean quadratic error (MQE):

and mean relative quadratic error (MRQE):

to evaluate the prediction ability of the model, where $\operatorname{Var}(Y)$ is the empirical variance of $\left\{Y_{t}\right\}_{t = 321}^{400}$. To demonstrate the superiority of the proposed model, we compare its prediction results with three standard models: the semi-functional linear model (SFL) from Zhou and Chen ^[10], the partial functional linear regression model (PFL) by Zhang et al. ^[26], and the partial functional linear regression model with autoregressive error (PFLAR) by Yu et al. ^[23]. Each simulation experiment is repeated 200 times, and the results are summarized in Table 3. From Table 3, we can observe that:

(ⅰ) When the true model is SFLAR, the proposed model has the best prediction accuracy.

(ⅱ) Using our proposed SFLAR to fit PFLAR data has the same performance as using PFLAR to fit PFLAR data.

5.   Applications to real data

In this section, the model and method proposed in this paper will be applied to study the power data set, which can be downloaded from https://www.eia.gov. This paper aims to forecast the monthly retail sales of electricity $C_i$ kilowatt hour (KWH) in the commercial sector. For this purpose, it is considered that the US monthly electricity consumed by commercial sectors which includes the power consumption data from January 1973 to January 2016 (517 months) (Figure 5) and their annual average retail price $P_t$ (per KWH, including tax) (43 years).

As shown in Figure 5, the time series shows some linear trends and some heterogeneity in the variance structure. Similar to Aneiros and Vieu ^[22] and Yu et al. ^[23], we eliminate the heteroscedasticity and the linear trend of the electricity retail sales data by differencing the $\ln$(data), and obtain the time series (see Figure 6): $\mathcal{X}_{i^{'}} = \ln \left(C_{i^{'}+1}\right)-\ln \left(C_{i^{'}}\right), i^{'} = 1, 2, \cdots, 516$. Let $X_t(s) = \left\{\mathcal{X}_{12(t-1)+s}, t = 1, 2, \cdots, 43, s = 1, 2, \cdots, 12 \right\}$ be the monthly $\ln$(difference electricity retail sales data).

There is only one observation per month, so we use ten cubic B-spline basis functions to convert the discrete monthly power data into functional continuous smooth annual curve predictors $\{X_t(u), u \in [1,12], t = 1, 2, \cdots, 43\}$ (see Figure 7).

The response variable is $Y_t(s) = \mathcal{X}_{12t + s}$, which represents the power consumption data in the $s$-month of the $t$-year, where $t = 1, 2, \cdots, 42, s = 1, 2, \cdots, 12$.

Similar to Aneiros and Vieu ^[22] and Yu et al. ^[23], the power data set is divided into the training sample $\left\{Y_{t}(s), P_{t}, X_{t}(u)\right\}_{t = 1}^{41}$ and the test sample $\left\{Y_{42}(s), P_{42}, X_{42}(u)\right\}$. According to Yang et al.^[24], we consider the following semi-functional linear model (SFLAR) with autoregressive error to fit the training set:

Then, the test set verifies the model's predicted value $Y_{42}(s)$. To demonstrate the superiority of the proposed model, we also compare its prediction results with three standard models: The semi-functional linear model (SFL), the partial functional linear regression model (PFL), and the partial functional linear regression model with autoregressive error (PFLAR). In addition, similar to Aneiros and Vieu ^[22] and Yu et al.^[23], we also utilize the mean quadratic error (MQE):

and mean relative quadratic error (MRQE):

to evaluate the prediction ability of the model, where $\operatorname{Var}(Y(s))$ is the empirical variance of $\left\{Y_{i}(s)\right\}_{i = 1}^{41}$.

The prediction effects of the four models under the test set are summarized in Table 4 and Figure 8. It is observed from Figure 8 that the four models exhibit consistency with the changing trend of the curve in Figure 7. Additionally, the prediction results from the test set in Table 4 and Figure 8 indicate that the SFL is closer to the actual value than the PFL. The MQE and MRQE of the SFL are also smaller than those of the PFL, suggesting that the prediction effect of the SFL is superior to the PFL, regardless of whether the autoregressive error is considered. Furthermore, the model with autoregressive error demonstrates a better prediction effect than the model without, implying that the former is more reasonable. The SFLAR model exhibits the smallest MQE and MRQE among the four models, signifying that the SFLAR model proposed in this paper yields the best prediction effect. Therefore, using the SFLAR model for this data is more reasonable.

6.   Discussion

The results of our study on forecasting monthly retail sales of electricity using an SFLAR model provide important insights into the dynamics of electricity consumption, which captures the influence of both functional and non-functional predictors of electricity sales. Incorporating autoregressive errors allowed the model to effectively capture temporal dependencies within the monthly sales data, thereby enhancing forecasting accuracy. However, our model also has some limitations. For example, we only consider one-dimensional covariate $Z$ in our model. Although our model can be easily generalized to multidimensional situations, it may encounter the curse of dimensionality and computational burden. As one of the anonymous reviewers said, a single index model could provide much more efficient results. Hence, we can extend our model to the semi-functional single index model in future research. At the same time, deep neural networks have been extensively used for nonparametric regression. This also constitutes a topic for our future research.

7.   Conclusions

This study focused on the SFLAR model, incorporating a random error sequence with an $AR(p)$ structure. By utilizing Mercer's theorem, the linear component of the function is dimensionally reduced using the method of functional principal component basis expansion. Meanwhile, the nonparametric function is approximated using the B-spline function. The slope function $\beta(\cdot)$, nonparametric function $g(\cdot)$, and autoregressive coefficient $a$ are estimated using a two-step iterative algorithm. Theoretical properties of the estimated parameters are provided under specific regular conditions. Based on the simulation results, the proposed model and estimation method perform well even with a limited sample size. The real data analysis also indicates that the SFLAR model has the best prediction accuracy among these four models when autoregressive errors exist. Therefore, the SFLAR model proposed in this paper is more suitable for this power data.

Author contributions

Bin Yang: Investigation, funding acquisition, formal analysis, validation, software, data curation, writing-original draft, writing-review and editing; Min Chen: Conceptualization, methodology, supervision, project administration, validation; Jianjun Zhou: Conceptualization, funding acquisition, investigation, resources, validation, data curation, writing-original draft, writing-review and editing. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors are grateful to the two anonymous referees and the editor for their constructive remarks that greatly improved the manuscript. The research is partially supported by a grant from the National Key R & D Program of China, Grant/Award Number: 2022YFA1003701, National Nature Science Foundation of China (No.11861074, 12261051), Applied Basic Research Project of Yunnan Province (2019FB138), Key Project of Yunnan Provincial Basic Research (Grant No.202401AS070061), and the Yunnan Provincial Department of Education Science Research Fund (Grant No.2024J0087).

Conflict of interest

The authors declare no conflict of interest.