A Functional Data Approach for Continuous-Time Analysis Subject to Modeling Discrepancy under Infill Asymptotics

Abstract

Parametric continuous-time analysis often entails derivations of continuous-time models from predefined discrete formulations. However, undetermined convergence rates of frequency-dependent parameters can result in ill-defined continuous-time limits, leading to modeling discrepancy, which impairs the reliability of fitting and forecasting. To circumvent this issue, we propose a simple solution based on functional data analysis (FDA) and truncated Taylor series expansions. It is demonstrated through a simulation study that our proposed method is superior—compared with misspecified parametric methods—in fitting and forecasting continuous-time stochastic processes, while the parametric method slightly dominates under correct specification, with comparable forecast errors to the FDA-based method. Due to its generally consistent and more robust performance against possible misspecification, the proposed FDA-based method is recommended in the presence of modeling discrepancy. Further, we apply the proposed method to predict the future return of the S&P 500, utilizing observations extracted from a latent continuous-time process, and show the practical efficacy of our approach in accurately discerning the underlying dynamics.

1. Introduction

Much of finance and economics is about the study of dynamics over time, in which analysis using time-series data plays a vital part. Despite discrete observations in practice, many time-series data such as stock prices, interest rates, and GDP are essentially drawn from their continuous-time underlying processes, for which consistent estimation can be achieved from continuous-time analyses but not necessarily from their discrete-time counterparts. Thus, the former has become increasingly incorporated by modern time-series analysis (e.g., [1,2,3,4]).

In traditional parametric continuous-time analyses, the modeling routinely originates in a discrete-time setting and is then extended to a continuous-time formulation as the length of the time window shrinks towards zero [3]. However, such derivation often imposes strong restrictions while specifying the convergence rates of frequency-dependent (hereafter, f.d.) parameters, resulting in models that are unadaptable to real-world data and thus forecasting failure [4,5,6,7,8,9,10,11,12]. For instance, with different sets of conditions regarding the converging speeds of the parameters, many GARCH-like processes have various diffusion processes (e.g., [6,9,11,13]). Moreover, through a testing procedure involving discrimination between the classes of deterministic and stochastic volatility models, ref. [7] discovered that, unlike the traditional stochastic view of short-term interest rates, the Japanese short-term interest rate indeed follows a deterministic pattern. Consequently, there arises discrepancy in continuous-time modeling, as it is debatable which assumption is the correct one (e.g., [14], pp. 176–178, [15]), and mistaking the assumption will lead to a misspecified limit and thereby unreliable analysis results. In such circumstances, nonparametric approaches appear to be appealing tools that can adapt to the true limits based on their data-driven nature and bypass this discrepancy in continuous-time analysis [16].

There is a large body of literature demonstrating the effectiveness of nonparametric methods for accurate estimation and forecasting in cases where parametric assumptions are deemed to be inadequate (e.g., [17,18,19,20,21,22,23,24,25,26]). Many studies incorporate nonparametric methods to approximate the density of the states in the absence of a closed-form expression, so that they could use maximum likelihood estimation (MLE) for continuous-time diffusion models (e.g., [22,27,28]). However, they still presume the parametric format of the underlying processes. Another stream of literature uses kernel-based methods to estimate the processes that are nonanticipative smooth functions with unknown structures (e.g., [23]) and make forecasts through conditional density estimation (e.g., [19,20,29,30]). It has been shown, however, that traditional kernel estimators can become inconsistent as the sampling density grows despite the underlying processes becoming further revealed [18,31], and the infill asymptotics (i.e., the asymptotic properties achieved as the sample becomes increasingly dense [32]) require the careful consideration of the dependence among observations, which is substantial work [33].

The main contribution of our paper is to propose a fitting and forecasting approach from the viewpoint of functional data analysis (FDA) to accommodate f.d. data structures and make good use of high-frequency data in order to achieve robust infill asymptotics. Indeed, the FDA fitting method, by design, deals with a collection of subsets of finite-dimensional parameter spaces, which becomes richer and denser with an enlarging sample, and the target function can be consistently estimated by optimizing an empirical criterion (e.g., [34,35,36]). For such a feature, in our method we employ the FDA-based approach to fit the underlying process using local polynomial bases and obtain the forecasts by extending the movement of the process based on the boundary derivatives of the functional fitting. FDA is renowned for its ability to uncover the dynamics of unknown continuous processes without necessitating stringent assumptions about the data structure [26,37,38]. However, to the best of our knowledge, its application in tracing the f.d. data structures and performing forecasting in continuous-time analysis remains underexplored in the existing literature. Thus, with a simulation study, this paper illustrates the procedure and properties of the proposed FDA-based method and emphasizes that this method excels when the convergence rate of f.d. parameters is undefined, effectively mitigating the common problem of misspecification prevalent in parametric estimations. In addition, we remark that functional approaches can be further developed in this context to achieve desired infill asymptotics with bounded or unbounded domains. To improve readability, all proofs are placed in Appendix A.

In addition to the simulation analysis, we substantiate the practical utility of the proposed method in forecasting S&P 500 prices, a topic of enduring significance that has consistently commanded the interest of financial practitioners and researchers within the field. Throughout the years, various statistical and mathematical methods have emerged to offer insights into stock price trends including but not limited to time-series analysis (e.g., [39,40,41]); machine learning algorithms (e.g., [22,42,43]); technical indicators (e.g., [44,45,46]); and fundamental analysis (e.g., [47,48,49]). However, to the best of our knowledge, the utilization of FDA for future stock price predictions appears to be infrequent. This paper adopts the FDA-based approach to predict the S&P 500 log-returns. Through the root mean squared forecast error (RMSFE) among different sample sizes and frequencies, we examine the performance of the proposed method in forecasting the process in a given step ahead.

The rest of this paper is structured as follows. In Section 2, we explain the FDA-based method, followed by a discussion of the large-sample properties. Section 3 sets up a simulation study and takes strong GARCH(1,1) as a motivating example to show that in continuous-time analysis, while parametric methods—such as MLE—encounter discrepancy, the proposed FDA-based method can provide robust and reliable estimation and forecasting. We then present an application of the proposed method in predicting future S&P 500 prices in Section 4, and Section 5 concludes the paper.

2. Methods

In this section, we explain the functional fitting and forecasting methods, followed by a discussion of the large-sample properties of the estimators. The method concerns target functions that are continuous but not necessarily differentiable. While smoothness assumptions are often imposed in FDA to accommodate incomplete and noisy observations, as noted by [50], we construct the functional estimator through the Bernstein polynomial approximation, which induces the convergence that requires only the continuity of the underlying process.

2.1. Fitting and Forecasting

Consider a complete probability space $(\mathsf{\Omega },\mathcal{F},P)$ and a re-scaled bounded time interval $[0,1]$. The sample path of the stochastic process $\left\{X\right(t):t\in [0,1\left]\right\}$ for state $\omega \in \mathsf{\Omega }$ over $[0,1]$ is denoted as $X_{\omega }$. Realizations are drawn at countably many time points with observational error ${ϵ}_{\omega ,t}$; as such,

$$\begin{array}{c}\hfill X_{\omega ,t}=X_{\omega }\left(t\right)+{ϵ}_{\omega ,t}.\end{array}$$

(1)

In parametric modeling, $X_{\omega }\left(t\right)$ is the conditional mean with a specified functional form, whereas in FDA, $X_{\omega }\left(t\right)$ is the underlying sample path assumed to be continuous at any given point t. Thus, in the discussion below, we establish the properties of our fitting and forecasting quantities through an arbitrarily close approximation of $X_{\omega }\left(t\right)$ by the Bernstein polynomial $B_K(t,X_{\omega })$, so that for each $K\in \mathbb{N}$,

$$\begin{array}{c}\hfill B_K(t,X_{\omega }):=\sum _{k=0}^K{X}_{\omega }\left(\frac{k}{K}\right)\times \left(\genfrac{}{}{0pt}{}{K}{k}\right)\times t^k{(1-t)}^{K-k},\end{array}$$

(2)

whence we have the following lemma.

Lemma 1.

Consider $B_K(t,X_{\omega })$ as in Equation (2) with the continuous sample path $X_{\omega }$ for any given $\omega \in \mathsf{\Omega }$. For any $\epsilon >0$, there exists some $\delta \left(\epsilon \right)>0$ that induces $|t_1-t_2|\le \delta \left(\epsilon \right)\Rightarrow |X_{\omega }\left(t_1\right)-X_{\omega }\left(t_2\right)|\le \epsilon /2$ for all $t_1,t_2\in [0,1]$. Then, for all $K\ge 4{sup}_{t\in [0,1]}\left|X_{\omega }\left(t\right)\right|/\left({\delta }^2\left(\epsilon \right)\epsilon \right)$, the following holds:

$$\underset{t\in [0,1]}{sup}|B_K(t,X_{\omega })-X_{\omega }\left(t\right)|\le \epsilon .$$

Lemma 1 suggests that, fixing everything else, for any given $\epsilon$, a high-degree polynomial is needed for approximation if the continuity of $X_{\omega }$ requires a very small $\delta \left(\epsilon \right)$, while a low-degree polynomial can achieve the desired approximation if a large $\delta \left(\epsilon \right)$ suffices for the given $\epsilon$ (see, e.g., Theorem 5.14 [51]). This lemma also implies that while fitting any continuous process, the typical smoothness assumptions for functional data can be imposed on the Bernstein approximating polynomials $B_K(t,X_{\omega })$, and the size of the remaining estimation error can be regulated by a high-degree $B_K(t,X_{\omega })$. There are three remarks following Lemma 1. First, this lemma is based on the continuity of $X_{\omega }\left(t\right)$, which does not require the continuity of the observed series. Second, the lemma holds for any continuous sample path and thus can also be applicable when $X_{\omega }\left(t\right)$ itself is a diffusion process of Hölder continuity. Last but not the least, there are various methods to obtain approximations to a continuous function with desirably small errors (Chapter 10 [52]), and here we adopt the Bernstein polynomial for the delivery of the discussion. Hence, the establishment of the asymptotics in the current study is achieved with the help of the parametrization of $B_K(t,X_{\omega })$, while in practice, we do not directly estimate or determine the degree K in $B_K(t,X_{\omega })$. Instead, the degree of the approximating polynomials will be controlled through the construction of the functional estimator, as explained below.

In the spirit of FDA [38], one can construct a B-spline representation of a certain degree for the polynomial approximation $B_K(t,X_{\omega })$, such that $B_K(t,X_{\omega })={\mathsf{\Phi }}^{⊺}(·)\mathit{C}$, where $\mathsf{\Phi }$ is a B-spline basis that contains a vector of Q basis functions, ${\mathit{C}}_{\omega }$ is a vector of Q corresponding basis coefficients, and the superscript “⊺” indicates the transpose operation. Then, Equation (1) can be rewritten as

$$\begin{array}{c}\hfill X_{\omega ,t}\approx {\mathsf{\Phi }}^{⊺}\left(t\right){\mathit{C}}_{\omega }+{ϵ}_{\omega ,t}.\end{array}$$

(3)

With the predetermined basis functions $\mathsf{\Phi }$ and the J observations ${\left\{X_{t_j}\right\}}_{j=1}^J:={\left\{X_{\omega ,t_j}\right\}}_{j=1}^J$ at ${\left\{t_j\right\}}_{j=1}^J\subseteq [0,1-\Delta ]$ for some $0<\Delta <1$, the coefficients ${\mathit{C}}_{\omega }$ are estimated by minimizing the penalized squares $J^{-1}{\sum }_{j=1}^J{[X_{t_j}-\tilde{X}\left(t_j\right)]}^2+\lambda {\int }_0^1{\tilde{X}}^{\left(2\right)}\left(t\right)dt$ with a tuning parameter $\lambda$ and some fitted function $\tilde{X}(·)$, which yields

$$\begin{array}{c}\hfill {\tilde{\mathit{C}}}_{\omega }:={\left[\frac{1}{J}\sum _{j=1}^J\mathsf{\Phi }\left(t_j\right){\mathsf{\Phi }}^{⊺}\left(t_j\right)+\lambda {\int }_0^1{\mathsf{\Phi }}^{\left(2\right)}\left(t\right){\left\{{\mathsf{\Phi }}^{\left(2\right)}\left(t\right)\right\}}^{⊺}dt\right]}^{-1}\frac{1}{J}\sum _{j=1}^J\mathsf{\Phi }\left(t_j\right)X_{\omega ,t_j},\end{array}$$

(4)

where the superscript “$\left(2\right)$” indicates the second-order derivative of a function (the second-order derivative of the fitted function is one option for the roughness penalty—other penalties can also be applied depending on the situation). In practice, Q and $\lambda$ can be selected through certain data-driven algorithms (e.g., [38]). Henceforth, the fitted process of $X_{\omega }\left(t\right)$ is defined as

$$\begin{array}{c}\hfill {\tilde{X}}_{\omega }\left(t\right):={\mathsf{\Phi }}^{⊺}\left(t\right){\tilde{\mathit{C}}}_{\omega },\forall t.\end{array}$$

(5)

By choosing a B-spline basis of a proper degree larger than some value R for the fitted process ${\tilde{X}}_{\omega }\left(t\right)$, the $\Delta$-step-ahead predictor at any given t, ${\widehat{X}}_{\omega }(t+\Delta )$ is defined as a truncated Taylor expansion with the first R derivatives of the fitted functions, such that

$$\begin{array}{c}\hfill {\widehat{X}}_{\omega }(t+\Delta ):=\sum _{r=0}^R\frac{1}{r!}{\Delta }^r{\tilde{X}}_{\omega }^{\left(r\right)}\left(t\right)=\sum _{r=0}^R\frac{1}{r!}{\Delta }^r{\tilde{\mathit{C}}}_{\omega }^{⊺}{\mathsf{\Phi }}^{\left(r\right)}\left(t\right)\forall t\in (0,1-\Delta ].\end{array}$$

(6)

The validity of the above Taylor expansion requires the fitted function ${\tilde{X}}_{\omega }\left(t\right)$ to be at least R times continuously differentiable, while the underlying process $X_{\omega }\left(t\right)$ only needs to be continuous. We want the relationship between the order R of continuous differentiability and the degree K of the approximated function $B_K(t,X_{\omega })$ to be such that $R=K+1$, as the Taylor polynomial approximation will have a non-zero residual if $B_K(t,X_{\omega })$ is of a higher degree (i.e., if $R-1<K$), while a higher degree of the functional estimator than necessary (i.e., $R-1>K$) will bring extra volatility to the fitting and forecasting values. In practice, R can also be selected through certain data-driven methods, such as that explained in the simulation and empirical studies below. Due to the features of the B-spline basis employed, one can adopt the order of $R+3$ to ensure the desired differentiability of ${\tilde{X}}_{\omega }\left(t\right)$. The half-open interval $(0,1-\Delta ]$ indicates that, technically, the forecast cannot be achieved without any observation, and the forecast can be obtained as far as $t+\Delta =1.$ Furthermore, although the derivatives ${\tilde{X}}_{\omega }^{\left(r\right)}\left(t\right)$ are applied in the approximation, we do not actually need the convergence of these derivatives to, even if they exist, the derivatives of the truth—only the convergence of the forecasting function in the level is of interest.

2.2. Large-Sample Properties

We now discuss the consistency and the asymptotic normality of our functional predictors based on the following assumptions.

Assumption 1.

Consider a continuous sample path $X_{\omega }\left(t\right)$ where for any $\epsilon >0$, there exists $\delta \left(\epsilon \right)>0$ such that $|t_1-t_2|\le \delta \left(\epsilon \right)$ implies $|X_{\omega }\left(t_1\right)-X_{\omega }\left(t_2\right)|\le \epsilon /2$. Then, we have the following:

(a) 

The sample path $X_{\omega }\left(t\right)$ is observed on a set of evenly spaced time points.

(b) 

The observational error ${ϵ}_{\omega ,t_j}$ is uncorrelated across ${\left\{t_j\right\}}_{j=1}^J\subseteq [0,1]$, with $\mathbb{E}[{ϵ}_{\omega ,t}\mid X_{\omega }\left(s\right)]=0$ for $s<t$ and $\mathit{Var}\left[{ϵ}_{\omega ,t}\right]<c<\infty$ for all $t\in [0,1]$ and some constant c.

(c) 

$Q\sim J^{{\alpha }_1}$ and $\lambda \sim J^{{\alpha }_2}$ with $0<{\alpha }_1<1$ and ${\alpha }_2<0$.

For simplicity, we assume in Assumption 1(a) that the observations are equally spaced without a loss of generality. A less restrictive version would be to assume that the ratio between the lengths of the largest and the smallest time windows is bounded above and away from zero as in [53]. Part (b) states a sufficient condition to achieve the consistency of the functional estimator ${\tilde{X}}_{\omega }$ as in Equation (5), while a proper “low correlation” assumption for this error term is also sufficient. Finally, part (c) indicates that the dimension of the basis expansion shall increase with the sample size to allow for a consistent functional estimator ${\tilde{X}}_{\omega }$, and, meanwhile, the tuning parameter shall be kept to $o\left(1\right)$ so that the estimation error introduced by the roughness penalty dies down as the sample size grows towards infinity.

For any given $0<\Delta <1$, we consider the convergence of the forecasting process ${\widehat{X}}_{\omega }$ over $(\Delta ,1)$ or, equivalently, over $t\in (0,1-\Delta )$. Note that even though the forecasting process can be defined at the point $t+\Delta =1$ as in Equation (6), the asymptotics exclude the point of $t=1-\Delta$ due to the limited boundary performance.

Theorem 1

(Convergence of ${\widehat{X}}_{\omega }$ given $\Delta$). Consider a continuous sample path $X_{\omega }$ on $[0,1]$ such that for any $\epsilon >0$, there exists $\delta \left(\epsilon \right)>0$, which induces $|t_1-t_2|\le \delta \left(\epsilon \right)\Rightarrow$$|X_{\omega }\left(t_1\right)-X_{\omega }\left(t_2\right)|\le \epsilon /2$.

Let ${\mathcal{F}}_K$ be the collection of all the (local) polynomials over $[0,1-\Delta ]$ that are at least $K+1$ times continuously differentiable, and $F_K:{\mathcal{F}}_K\to \mathbb{R}$ be the linear functional such that $F_K\left(\psi \right):=〈{\widehat{X}}_{\omega }(·+\Delta )-B_K(·+\Delta ,X_{\omega }),\psi 〉$ for $\psi \in {\mathcal{F}}_K$.

If Assumption 1 holds, then for any given K and $R=K+1$, $F_K\left(\psi \right)\stackrel{p}{⟶}0$ as $J\to \infty$ for all $\psi \in {\mathcal{F}}_K$, whence under $K\to \infty$, ${sup}_{t\in (0,1-\Delta )}|{\widehat{X}}_{\omega }(t+\Delta )-X_{\omega }(t+\Delta )|\stackrel{p}{⟶}0$ as $J\to \infty$.

Theorem 1 establishes the convergence of the forecasting process ${\widehat{X}}_{\omega }(·+\Delta )$ through its weak convergence to the Bernstein approximating polynomial $B_K(·+\Delta ,X_{\omega })$ and the uniform convergence of $B_K(·,X_{\omega })$ to $X_{\omega }$. Specifically, for any given K and $R=K+1$, ${\widehat{X}}_{\omega }(·+\Delta )$ weakly converges to the approximating polynomial $B_K(·+\Delta ,X_{\omega })$ with the enlarging sample size J, while $B_K(·,X_{\omega })$ uniformly converges to the true path $X_{\omega }$ with an enlarging K and regardless of the sample size J. Henceforth, ${\widehat{X}}_{\omega }(·+\Delta )$ converges uniformly to the true path $X_{\omega }$ with an enlarging K and J. For $R=K+1$ with $K\ge 4{sup}_{t\in [0,1]}\left|X_{\omega }\left(t\right)\right|/\left({\delta }^2\left(\epsilon \right)\epsilon \right)$, it is worth noting that for a B-spline basis of order $R+3$ to span the function space, the number of basis functions Q needs to be such that $Q\ge 2(R+2)+1$, while Q grows towards infinity at a lower speed than J; hence, we also have $K=o\left(J\right)$ as $K,J\to \infty$.

Now we explore the asymptotic distribution of the functional prediction ${\widehat{X}}_{\omega }$. However, before we state any further assumptions or the resulting asymptotic properties, it is important to note that the realizations are subject to observational error ${ϵ}_{\omega ,t}$, the behavior of which can impact the asymptotic distribution of ${\widehat{X}}_{\omega }$. For example, when a continuous sample path $X_{\omega }$ is observed without noise, it is more plausible to consider the functional fitting as an interpolation—in this case, we do not necessarily have asymptotic normality for ${\widehat{X}}_{\omega }$. On the other hand, if $X_{\omega }$ is observed with noise, we can establish pointwise asymptotic normality under the following sufficient conditions:

Assumption 2.

(a) 

The noise ${ϵ}_{t_j}$ is independent across all ${\left\{t_j\right\}}_{j=1}^J\subseteq [0,1]$.

(b) 

$Q\sim J^{{\alpha }_1}$ and $\lambda \sim J^{{\alpha }_2}$ with ${\alpha }_1>0$ and ${\alpha }_2<0$, such that $Q^2\lambda =o\left(J^{-1/2}\right)$ and $J^{-1}{Q}^2=O\left(J^{-1/2}\right)$.

Assumption 2(a) states a sufficient condition for generating the asymptotic normality by the Lyapunov CLT, which allows for non-identically distributed variables. If one has mixing random processes, a different version of this assumption can be specified by imposing identical distribution, where a different CLT can be applied without affecting the result regarding asymptotic normality. Part (b) assumes a stronger condition on the orders of the two parameters Q and $\lambda$, so that the non-normal part of the estimation error is dominated and will not be inflated while deriving the asymptotic distribution.

Theorem 2.

Consider a continuous sample path $X_{\omega }$ on $[0,1]$ where for any $\epsilon >0$, there exists $\delta \left(\epsilon \right)>0$ such that $|t_1-t_2|\le \delta \left(\epsilon \right)$ implies $|X_{\omega }\left(t_1\right)-X_{\omega }\left(t_2\right)|\le \epsilon /2$.

Under Assumptions 1 and 2, the following asymptotic normality holds for all $t\in (0,1-\Delta ]$, $R\to \infty$, and $J\to \infty$:

$$\begin{array}{c}\hfill {\left[V(t;\Delta ,R,\mathsf{\Phi })\right]}^{-1/2}\left[{\widehat{X}}_{\omega }(t+\Delta )-X_{\omega }(t+\Delta )\right]\stackrel{d}{⟶}\mathcal{N}(0,1),\end{array}$$

where ${\sigma }_j^2:=\mathit{Var}\left[{ϵ}_{t_j}\right]$ for all j, and

$$\begin{array}{c}V(t;\Delta ,R,\mathsf{\Phi }):=\frac{1}{J^2}\sum _{j=1}^J{\sigma }_j^2{A}_{t_j}^2(t;\Delta ,R,\mathsf{\Phi }),\hfill \\ A_{t_j}(t;\Delta ,R,\mathsf{\Phi }):=\sum _{r=0}^R\frac{1}{r!}{\Delta }^r{\left\{{\mathsf{\Phi }}^{\left(r\right)}\left(t\right)\right\}}^{⊺}{\left\{{\int }_0^1\mathsf{\Phi }\left(s\right){\mathsf{\Phi }}^{⊺}\left(s\right)ds\right\}}^{-1}\mathsf{\Phi }\left(t_j\right).\hfill \end{array}$$

Theorem 2 constructs the pointwise asymptotic normality of the functional predictor ${\widehat{X}}_{\omega }$ through the asymptotically normal prediction error from ${\widehat{X}}_{\omega }$ to the Bernstein approximating polynomial $B_K(t,X_{\omega })$ and the convergence of $B_K(t,X_{\omega })$ to $X_{\omega }$ as K increases towards infinity. Specifically, ${\widehat{X}}_{\omega }(t+\Delta )-X_{\omega }(t+\Delta )$ can be decomposed into ${\widehat{X}}_{\omega }(t+\Delta )-B_K(t+\Delta ,X_{\omega })$ and $B_K(t+\Delta ,X_{\omega })-X_{\omega }(t+\Delta )$, where the former is inflated by $V(t;\Delta ,R,\mathsf{\Phi })$ to an asymptotic normal, and the latter is dominated because a desirably small approximation error $B_K(t+\Delta ,X_{\omega })-X_{\omega }(t+\Delta )$ can be achieved by a sufficiently large K (and R). Meanwhile, since the inverse factorial series ${\sum }_{r=0}^{\infty }{(r!)}^{-1}$ converges, when $R\to \infty$ we have ${\sum }_{r=0}^R{(r!)}^{-1}{\Delta }^r<C<\infty$ for some small $\Delta$ and some fixed $C\in \mathbb{R}$, whence $A_{t_j}(t;\Delta ,R,\mathsf{\Phi })$, and thus $A_{t_j}^2(t;\Delta ,R,\mathsf{\Phi })$ is bounded. Then, by Assumptions 1(b) and 2(a), the term $V(t;\Delta ,R,\mathsf{\Phi })$—and so is the variance ${\widehat{X}}_{\omega }(t+\Delta )-X_{\omega }(t+\Delta )$—is of order $J^{-1}$. As a result, a $\sqrt{J}$-asymptotic normality can be achieved.

3. Simulation

We now explain the numerical procedure of the proposed methods through a simulation study and compare its performance to that of a parametric approach under both correct specification and misspecification, thereby revealing FDA’s superiority in tracking the f.d. data structures of stochastic processes. The simulation study was conducted using the software R (version 4.0.5) and MATLAB (version R2023a), and the computer code is available in the Supplementary Materials.

3.1. The Data-Generating Process

The data-generating process is based on the strong GARCH (1,1), where $y_k=(S_k-S_{k-1})/S_{k-1}$, $k=1,2,\dots$ is the arithmetic return on a financial asset with the price $S_k$. Let h be the time window, and recall that a strong GARCH (1,1) process is represented by $y_k=\mu +{ϵ}_k$ with ${ϵ}_k\sim \mathcal{N}(0,V_k)$ given the $\sigma$-algebra generated by ${ϵ}_{k-1}$ and $V_k={\omega }_h+{\xi }_h{ϵ}_{k-1}^2+{\gamma }_h{V}_{k-1}$ for the f.d. parameters ${\omega }_h,{\xi }_h,{\gamma }_h>0$ and ${\xi }_h+{\gamma }_h<1$ [54,55].

As h shrinks, $\left(S_{kh},V_{kh}\right)$ weakly converges to its continuous-time limit $\left(S_t,V_t\right)$ [6]. This is a simple but effective example, since given different assumptions as to the f.d. parameters’ convergence rates (i.e., the convergence of ${\xi }_h$ can have rate $\sqrt{h}$ or rate h), $\left(S_t,V_t\right)$ has been shown to be a diffusion process solution to either a stochastic volatility (SV) model [13,56]

$$\begin{array}{c}\hfill \left[\begin{array}{c}d{s}_t\\ d{v}_t\end{array}\right]=\left[\begin{array}{c}a\\ \left(\beta -\frac{1}{2}{\sigma }^2\right)+\alpha exp(-v_t)\end{array}\right]dt+\left[\begin{array}{cc}\sqrt{1-{\rho }^2}exp\left(\frac{v_t}{2}\right)& \rho exp\left(\frac{v_t}{2}\right)\\ 0& \sigma \end{array}\right]\left[\begin{array}{c}d{W}_t^1\\ d{W}_t^2\end{array}\right],\end{array}$$

(7)

where $s_t:=log\left(S_t\right)$, $v_t:=log\left(V_t\right)$, and $W_t^1$ and $W_t^2$ are independent standard Brownian motions, or a deterministic volatility (DV) model [6]

$$\begin{array}{c}\hfill \left[\begin{array}{c}d{s}_t\\ d{v}_t\end{array}\right]=\left[\begin{array}{c}a\\ \beta +\alpha exp(-v_t)\end{array}\right]dt+\left[\begin{array}{c}exp\left(\frac{v_t}{2}\right)\\ 0\end{array}\right]d{W}_t^1.\end{array}$$

(8)

Indeed, despite the fact that the ARCH-type diffusion models under both SV and DV have been frequently applied in the literature to estimate continuous-time processes (e.g., [8,9,22,27,57,58]), there has been considerable debate as to the choice of the convergence rate assumptions, which makes parametric analysis such as MLE tailored for either limit questionable (e.g., [14,59]). More generally, if the f.d. data structures have different limiting processes under different convergence conditions, it would be impractical for parametric continuous-time generalizations to exhaust all possible limits to choose the correct one.

To demonstrate how the FDA-based method can bypass this issue, we generate pseudo-continuous return and volatility processes, as shown in Equations (7) and (8), over the interval of $\left[0,1\right]$, resembling the “continuous-time” processes of five equally spaced trading points per day for one year of 252 trading days (i.e., 5 time points per day × 252 days per year × 1 year = 1260 data points in total). The true values of the parameters in

Table 1. Estimating SV parameters with MLE.

Parameter	a	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\sigma }$	$\mathit{\rho }$
True value	0.1	0.2	−8.5	2.7	−0.8
SV process fitted by an SV model (correct specification)
Estimate	0.004	0.289	−9.546	2.699	−0.800
Bias	−0.096	0.089	−1.046	−0.001	0.000
Rejection rate	0.066	0.080	0.047	0.056	0.062
SV process fitted by a DV model (misspecification)
Estimate	0.101	232.294	−259.989	–	–
Bias	0.001	232.094	−251.489	–	–
Rejection rate	0.060	0.789	0.902	–	–

Note: a 1000-time simulation allows for a “$\pm \sqrt{0.95\ast 0.05/1000}\ast 2=0.014$” error on the $5\%$ rejection rate under correct specification.

and

Table 2. Estimating DV parameters with MLE.

Parameter	a	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\sigma }$	$\mathit{\rho }$
True value	0.1	0.2	−8.5	0	0
DV process fitted by a DV model (correct specification)
Estimate	0.107	0.200	−8.531	–	–
Bias	0.007	0.000	−0.031	–	–
Rejection rate	0.048	0.051	0.053	–	–
DV process fitted by an SV model (misspecification)
Estimate	1.033	0.200	−8.500	0.000	0.154
Bias	0.933	0.000	0.000	0.000	0.154
Rejection rate	0.346	1.000	1.000	1.000	0.877

Note: a 1000-time simulation allows for a “$\pm \sqrt{0.95\ast 0.05/1000}\ast 2=0.014$” error on the $5\%$ rejection rate under correct specification.

are used to generate $N=1000$ pseudo-continuous log-return and log-volatility trajectories from the SV limit, denoted respectively by $s_S\left(t\right)$ and $v_S\left(t\right)$, and those trajectories from the DV limit, denoted respectively by $s_D\left(t\right)$ and $v_D\left(t\right)$. It should be noted that the volatility is not always observed in practice, and discussions addressing parametric analysis with such unobservability have been provided in the literature (e.g., [22,27,60]). However, the FDA-based method can handle the return and volatility processes separately; hence, in what follows, we assume that volatility is accessible without harming the analysis of the return. Should we be interested in the unobserved volatility, the same FDA-based method could be applied to suitable proxies. We take the pseudo-continuous-time processes as the observations, denoted as $s_S\left(t_j\right)$, $v_S\left(t_j\right)$, $s_D\left(t_j\right)$, and $v_D\left(t_j\right)$ for $j=1,\dots ,J$, to numerically mimic the scenario where the sampling time window becomes arbitrarily small, and the estimation is performed based on one-month rolling windows (i.e., $J=1260/12=105$) and eight-month rolling windows (i.e., $J=1260/12\ast 8=840$), respectively. Cases with an evenly spaced “daily” sampling density are also included to demonstrate the performance of the FDA-based method under different sample sizes and frequencies.

3.2. Fitting and Forecasting with FDA

Now we present the procedures for fitting and forecasting with FDA in steps (a) to (c) below, and the forecast evaluation is discussed in steps (d) and (e). Note that only the implementations in terms of $s_S\left(t\right)$ and $s_D\left(t\right)$ are explained here, and similar procedures can be applied to the $v_S\left(t\right)$ and $v_D\left(t\right)$ processes.

(a)

Provide a candidate pool for the vector $(\lambda ,R)$ of the tuning parameter and the Taylor expansion order. For each of the 1000 replications, take the length $J-10$ rolling windows (except for the one-month daily sample, where we take the length $J-3$ rolling windows) from the first J observed data points and estimate the underlying processes $s_S\left(t\right)$ and $s_D\left(t\right)$ as in Equations (4) and (5) using the B-spline basis of order $R+3$ and the number of basis functions $min\left\{2\right(R+2)+1,J+R+3\}$. We then obtain ${\tilde{s}}_S\left(t\right)$ and ${\tilde{s}}_D\left(t\right)$, respectively, for each rolling window with each pair of $(\lambda ,R)$ candidates.

(b)

For the first J observed data points from each of the 1000 replicates, compute the forecasting values ${\widehat{s}}_S(t+\Delta )$ and ${\widehat{s}}_D(t+\Delta )$ for each of the length $J-10$ rolling windows with each pair of $(\lambda ,R)$ candidates using Equation (6). The $(\lambda ,R)$ pair that minimizes the RMSFE over the ten rolling windows across all 1000 replications is selected and denoted by $(\widehat{\lambda },\widehat{R})$ to be used for later fitting and forecasting.

(c)

For each of the 1000 replications, perform the fitting and forecasting on the length J rolling windows (either one-month or eight-month) with the given observation frequency using the selected parameters $(\widehat{\lambda },\widehat{R})$, the basis order of $R+3$, and the number of basis functions $min\left\{2\right(R+2)+1,J+R+3\}$, yielding the rolling-window forecasting values ${\widehat{s}}_S(t+\Delta )$ and ${\widehat{s}}_D(t+\Delta )$. It is worth noting that the forecasting step of this method can be flexible, as $\Delta$ is technically a continuous quantity. In this simulation study, we only consider the forecasting step that is the same as the sampling time window, which we later refer to as the “one-step-ahead forecast”, while in general, for example, one can make a one-hour-ahead forecast with daily data or a one-day-ahead forecast with hourly data.

(d)

Implement the Kolmogorov-Smirnov ($K\text{-}S$) test on the pairs “${\widehat{s}}_S(t+\Delta )$ and $S_S(t+\Delta )$”, “${\widehat{s}}_S(t+\Delta )$ and $S_D(t+\Delta )$”, “${\widehat{s}}_D(t+\Delta )$ and $S_S(t+\Delta )$”, and “${\widehat{s}}_D(t+\Delta )$ and $S_D(t+\Delta )$” for all time points t to check whether the FDA-based predictors can correctly distinguish the true underlying processes in terms of distributions.

(e)

Calculate the RMSFE according to [61,62] between the pairs “${\widehat{s}}_S(t+\Delta )$ and $s_S(t+\Delta )$”, “${\widehat{s}}_S(t+\Delta )$ and $s_D(t+\Delta )$”, “${\widehat{s}}_D(t+\Delta )$ and $s_S(t+\Delta )$”, and “${\widehat{s}}_D(t+\Delta )$ and $s_D(t+\Delta )$” for the out-of-sample performance evaluation.

In practice, the relevant parameters can be selected more exhaustively for each forecast and each replication as needed; however, developing the asymptotic properties of such data-driven methods is not the focus of the current paper, so we simplify the search in step (a) by having only three possible Rs and three possible $\lambda$s for each search, forming nine pairs of $(\lambda ,R)$ candidates, and the resulting $(\widehat{\lambda },\widehat{R})$ is used for all the replicates in all the later forecasts. Also, the upper bound $J+R+3$ of the number of basis functions Q is used to simplify the search by first fixing Q at the sample size J while extending the fitting function beyond the fitting domain by $R+3$ equally spaced basis functions to capture the smoothness on $[1-\Delta ,1]$ and to avoid the boundary noise induced by the last $R+3$ (i.e., the order of the basis in our setting) basis functions.

3.3. Comparison with Parametric Methods

We adopted MLE as a benchmark parametric method for a performance comparison with our FDA approach as it is a commonly used approach for GARCH-like models (see e.g., [22,63]). It is imperative to underscore that within this context, MLE was employed as the benchmark parametric approach for the estimation of a sequence of discrete observations drawing from an unknown underlying diffusion process. While alternative parametric methodologies, like the general method of moments (GMM), are available for application, it is noteworthy that they all necessitate certain parametric assumptions about the data distribution, which could potentially result in erroneous estimation results and forecasts when dealing with f.d. data structures under infill asymptotics. Based on the stochastic differential equations of the two limits in Equations (7) and (8), the likelihood functions were obtained utilizing the (joint) normality of $d{W}_t^1$ and $d{W}_t^2$, and the differentials $ds$ and $dv$ were approximated by the corresponding first difference of the discrete observations. In favor of the parametric methods, we took the eight-month pseudo-continuous rolling windows as the observed series, and the comparison between the two approaches was conducted in terms of RMSFEs.

Thus, based on Equations (7) and (8), the forecasts based on conditional expectation are such that

$$\left\{\begin{array}{c}{\widehat{s}}_{t+\Delta }:=\widehat{\mathbb{E}}\left[s_{t+\Delta }|s_t\right]=s_t+\widehat{a}\Delta ,\hfill \\ {\widehat{v}}_{t+\Delta }:=\widehat{\mathbb{E}}\left[v_{t+\Delta }|v_t\right]=v_t+\left[\left(\widehat{\beta }-\frac{1}{2}{\widehat{\sigma }}^2\right)+\widehat{\alpha }exp(-v_t)\right]\Delta \hfill \end{array}\right.$$

(9)

for the SV limit and

$$\left\{\begin{array}{c}{\widehat{s}}_{t+\Delta }:=\widehat{\mathbb{E}}\left[s_{t+\Delta }|s_t\right]=s_t+\widehat{a}\Delta ,\hfill \\ {\widehat{v}}_{t+\Delta }:=\widehat{\mathbb{E}}\left[v_{t+\Delta }|v_t\right]=v_t+\left[\widehat{\beta }+\widehat{\alpha }exp(-v_t)\right]\Delta \hfill \end{array}\right.$$

(10)

for the DV limit, where $\widehat{a},\widehat{\beta },\widehat{\alpha },$ and $\widehat{\sigma }$ are the estimated parameters from MLE, and $\Delta$ is the forecasting step as previously defined.

3.4. Results

Figure 1 and Figure 2 present the forecasting results using the FDA-based method, where the underlying processes are indicated by gray lines, the one-step-ahead rolling forecasts are indicated by black lines, and the selected order of Taylor expansion is noted in the labels of the plots. The figures reveal that the proposed functional approach accurately traced the movements of both the underlying return and volatility processes.

Then, we used the $K\text{-}S$ test to compare the predicted and the underlying processes. The p-values of the test on the null of distribution equality over the entire time domain are presented in Figure 3 and Figure 4, where the dashed lines indicate the $5\%$ significance level. As depicted in the two plots from the first column of Figure 3, the true return and volatility underlying processes generated from the SV and DV models exhibited significantly disparate distributions at almost every fixed observation time $t_j$. The second and third columns of Figure 3 suggest that the null hypothesis of distribution equality was rejected when the underlying processes showed significant differences in the cross-comparison, while the FDA forecast shared the same distribution as its true underlying process at almost every $t_j$, as evident in Figure 4. The results imply that the one-step-ahead rolling forecast using the FDA-based method could preserve the distributions of the true underlying processes pointwisely and could distinguish the true underlying processes from a falsely assumed continuous-time limit in terms of their distributions.

The proposed method was further applied to a lower observation frequency and a smaller sample size, where it exhibited decent performance with short and sparse observed series; the corresponding results are illustrated in Appendix B.1. A summary of comparisons across different sample sizes and frequencies in terms of RMSFE is presented in Figure 5, where “1 Md” denotes the one-month rolling window with daily observations, “8 Md” denotes the eight-month rolling window with daily observations, “1 Mps” denotes the one-month rolling window with pseudo-continuous observations, and “8 Mps” denotes the eight-month rolling window with pseudo-continuous observations. The small yet discernible differences in the RMSFEs among the different domain spans and different sampling densities indicates that the FDA-based method benefited from both the infill sample growth and the increasing-domain sample growth, though more from the former than from the latter.

The results for the first rolling window of the MLE are shown in Table 1 and Table 2 to provide a snapshot of the performance of MLE, while the results for all 420 rolling-window estimations are summarized by boxplots in Appendix B.2. Unsurprisingly, with the misspecified models, the estimates of the affected parameters mostly departed from the true values with statistical significance, while MLE showed good asymptotic performance under the correct specification. As noted in Table 1 and Table 2, the rejection rates for some parameters slightly exceeded the reasonable range of errors induced by a binomial distribution with $N=1000$ trials and a probability of success at $5\%$ even when the model was correctly specified, which could be attributed to the insufficient number of observations. Nevertheless, in general, methods that can circumvent such modeling discrepancy will be useful to avoid misleading estimation results in continuous-time analysis.

Finally, to compare the performance of the FDA and MLE methods, we applied the MLE rolling-window estimation to perform one-step-ahead out-of-sample forecasting and compared the distributions of the relative RMSFEs (rRMSFEs) calculated as the ratio of the RMSFEs between FDA and MLE—a ratio below one signified that the FDA method provided a more precise forecast than the MLE method. The results for cases in which the MLE forecast was affected by the misspecification are shown in Figure 6, and those for cases in which the MLE forecast was unaffected are summarized in Figure A13 in Appendix B.2.

First, note that misspecifying an SV process and estimating it with a DV model did not cause significant bias in the estimation of the parameter a, which is the only parameter needed for the forecast of the return process, according to Equations (9) and (10). Thus, the estimation of $s_S$ was not significantly affected by misspecifying it as $s_D$, and the distribution of the rRMSFEs is summarized in Figure A13. However, the estimations for the parameters $\alpha$ and $\beta$ were affected by such a misspecification, which in turn caused inaccuracy and instability in the forecast of the volatility process, as shown in the first plot of Figure 6. Specifically, 47.9% of the rRMSFEs were lower than one, and around 44% of the rRMSFEs were lower than 0.01, while when MLE performed better, the rRMSFEs of the two methods were still within a comparable range. Hence, the FDA-based method appeared to be more stable and generally reliable. On the contrary, when a DV process was misspecified and estimated by an SV model, the estimation of all the parameters was affected statistically, but for the parameters $\alpha$, $\beta$, and $\sigma$, the biases were negligibly small. Hence, according to Equations (9) and (10) with $\widehat{\sigma }\approx 0$, the forecast for the volatility process would not be substantially affected, as summarized in Figure A13. For the forecast of the return process that relied on $\widehat{a}$, the estimates significantly deviated from the true value. However, since the forecast error was determined by the re-scaled estimation error $(\widehat{a}-a)\delta$ and the noise, the sizes of which were of the same order numerically, the final forecast error was not notably influenced by the error in $\widehat{a}$, and the forecast exhibited a similar performance to the FDA-based method, as shown in the second plot of Figure 6. However, in practice, when the step $\Delta$ became larger, the error in $\widehat{a}$ began to show a further impact on the forecast. Finally, the FDA-based method could still outperform the parametric method even in a correctly specified case, as shown in the last plot of Figure 6.

Above all, the FDA-based method showed general consistency and robustness against modeling discrepancy, while the parametric methods could demonstrate unreliable performance due to misspecification, and even when the parametric method excelled under correct specification, the forecast errors of the two methods still lay within a comparable range, as shown in Figure 6 and Figure A13. When DV was misspecified as SV, the MLE forecast for the volatility process surpassed FDA the furthest in terms of the RMSFE. However, for the forecast of this non-random function, both methods performed well, as shown by the RMSFEs in Figure 5. Hence, the proposed FDA-based method is recommended in the presence modeling discrepancy.

4. Empirical Study

This section illustrates the application of our FDA-based method in stock price prediction, through which we obtained the closing prices of the S&P 500 Index from the Bloomberg terminal. Predicting volatility entails similar processes to forecasting future price movements, for which historical volatility indices like the VIX could be incorporated into our methodology. In accordance with the procedure utilized in our prior simulation study, we proceeded to assess the predictive efficacy of the FDA method utilizing observations spanning a duration of one year. Specifically, one-month and eight-month rolling windows of hourly data (every hour from 9:30 to 15:30) and daily data (every day at 12:30) were employed, and the fitting and forecasting with parameter selection followed the procedures described in steps (a) to (c) from Section 3.

Figure 7 and Figure 8 depict the functional predictions across various sampling frequencies and rolling windows. Figure 9 and Figure 10 present the distribution of the ratio of absolute forecast errors to log-returns in the four corresponding scenarios. The percentage of the forecast errors was consistently small, suggesting a close alignment between the forecast and the actual log-returns, which underscored the robustness and reliability of the proposed FDA-based method.

5. Conclusions

In continuous-time modeling, parametric methods fail to provide reliable analysis when there is discrepancy due to the existence of multiple limits. This paper adopted FDA to uncover the true continuous-time underlying processes subject to f.d. data structures under infill asymptotics and suggested a forecasting method by integrating FDA with Taylor series expansion, also exploring an application of FDA in out-of-sample prediction.

Our theorems demonstrate that the FDA-based method only requires the continuity of the underlying sample path and the low correlation of the observation errors, and with proper basis expansions, enlarging samples with an increasing sampling density ensure that the functional estimator converges to a unique and well-defined limit. The simulation analysis showed that the FDA-based method is capable of distinguishing processes with different continuous-time limits in out-of-sample prediction. Furthermore, although the parametric method dominated under correct specification, the forecast errors of the two methods were still within a comparable range, while the FDA-based method showed a generally consistent and much more robust performance against possible misspecification when the parametric method was subject to modeling discrepancy, which makes the functional method a preferable tool for fitting and forecasting when there is uncertainty in modeling the underlying process. Further, we validated the practical applicability of the proposed method in predicting future S&P 500 prices. Our findings demonstrated that the functional approach yielded dependable predictions even in cases where the underlying process was unknown. It is important to note that, as shown in the empirical study, the suggested FDA-based method is not limited to any particular continuous-time model; instead, it is able to uncover the underlying process of many discretely observed stochastic processes, including equity returns.

Lastly, the proposed FDA-based method shows limitations in certain circumstances. First, as the most natural way to obtain the optimal values of the parameters determining the functional estimators and forecasts, e.g., the number and the order of the basis function, the tuning parameter $\lambda$, and the order of the Taylor expansion, is by data-driven algorithms, the search within such an optimization process can induce computationally expensive tasks. Additionally, in the presence of a continuous but volatile underlying process, our approach also necessitates increased computational resources for determining the high-order derivatives of the fitted functions. Moreover, the nature of the method implies that, with everything else fixed, it generally works more comfortably and easily for smoother processes than volatile ones and is more suited for capturing long-term patterns than short-term shocks. Therefore, the performance of the method may heavily rely on high-frequency data while fitting and forecasting highly volatile processes. For future research, more efficient methods to obtain the optimal values of the parameters are to be developed, and the effects of the estimated parameters on the overall performance of the fitting and forecasting are to be explored. Furthermore, the method’s performance with volatile underlying processes and/or sparse observations is to be further studied and improved.