An operator-level ARCH Model

An operator-level ARCH Model

\affil[1]{Department of Statistics, University of California, Davis, CA 95616 Davis, USA, \textsuperscript{\ddag}[email removed]} \affil[2]{Department of Mathematics, Ruhr University Bochum, 44780 Bochum, Germany, \textsuperscript{*}[email removed]} \affil[3,4]{Department of Statistics and Actuarial Science, University of Waterloo, ON N2L 0A4 Waterloo, Canada, \textsuperscript{\dag}[email removed], [email removed]}

{ MSC 2020 subject classifications: 60G10, 62F12, 62R10}

{ Keywords: ARCH; financial time series; functional data; parameter estimation; stationary solutions}

Introduction

Conditionally heteroscedastic processes are fundamental in financial modeling. To capture their dynamics, Engle1982 introduced the Autoregressive Conditional Heteroscedastic (ARCH) model, later extended by Bollerslev1986 to the Generalized ARCH (GARCH) model. Extensive reviews of uni- and multivariate (G)ARCH as well as other related volatility models can be found in AndersenEtAl2009, mvgarch2, mvgarch1, FrancqZakoian2019, and Gourieroux1997.

With advances in data processing and the growing availability of high-frequency financial data, functional volatility and (G)ARCH (f(G)ARCH) models have gained prominence. As a motivating example, consider observing on consecutive days $k\in\{1, \dots,N\}$ the price of the S&P 500 index at a resolution of one minute. Conveniently, these observations can be viewed as a sample of curves or functions $\{P_k(t)$, $i \in \{1,\dots,N\},\;t\in[0,1]\}$, where intraday time is normalized to the unit interval; see the left panel of Figure (ref). A transformation of particular interest are the overnight intraday log-returns $X_k(t)= \log P_k(t) - \log P_{k-1}(1)$; see the right panel of Figure (ref). Representation of the data as functions constructed from noisy high-dimensional intraday observations allows for the use of techniques from functional data analysis. We refer the reader to HorvathKokoszka2012 and HsingEubank2015 for introductions to functional data analysis, and to Bosq2000 and BosqBlanke2007 for introductions to linear time series processes in function spaces.

In seminal work, Hoermannetal2013 introduced the fARCH$(1)$ process to model conditional heteroscedasticity in functional time series data. Their model was extended to fGARCH$(1,1)$ by Aueetal2017, and to general fGARCH$(p,q)$ by Ceroveckietal2019. For further developments and modified functional GARCH and volatility models, see Andersen02042024, KearneyShangZhao2023, Kuehnert2019, Kuehnert2020, Laksacietal2025, LiSunLiu2025, RiceEtAL2023, SunYu2020.

Henceforth, we refer to these existing models as “pointwise” f(G)ARCH (pw-f(G)ARCH) models, since they are formulated as follows: For the function spaces $C[0,1]$ and $L^2[0,1]$ of continuous and square-integrable functions on $[0,1],$ the pw-fARCH$(p)$ model as introduced in Ceroveckietal2019 takes the form

Here, the model parameters are $\delta$, a positive function, and $\alpha_i$, bounded linear operators mapping non-negative functions to non-negative functions, with the $\varepsilon_k$ being i.i.d. and centered innovations. When there exists a stationary and causal solution to (ref) with finite second-order moments, the condition $\operatorname{\mathds E}(\varepsilon^2_0(t))=1$ identifies $\sigma_k^2(t)= \mbox{Var}(X_k(t)| \mathcal{F}_{k-1})$ as the pointwise conditional variance, where $\mathcal{F}_{k} =\sigma( \varepsilon_i, \; i \le k)$ is the information filtration generated by the innovations.

Despite its success in numerous applications, this model has some limitations. First, it only explicitly models the serial dependence structure of the pointwise conditional variance $\mbox{Var}(X_k(t)| \mathcal{F}_{k-1})$. For many tasks, for example the construction of a prediction set for $X_k$, one wishes to estimate the full conditional covariance kernel/operator of the sequence. Under model (ref), so long as it is well-defined, the conditional covariance kernel takes the form $$ \mbox{Cov}\big(X_k(t),X_k(s) \big| \mathcal{F}_{k-1}\big) \;=\; \operatorname{\mathds E}\!\big[ X_k(t) X_k(s) \big| \mathcal{F}_{k-1}\big] \;=\; \sigma_k(t) \sigma_k(s)C_\varepsilon(t,s), $$ where $C_\varepsilon$ is the covariance kernel of the innovations. As such, the conditional covariance forecast implied by model (ref) constitutes a “rank-one update” to the covariance of the innovations, which may be overly simplistic.

Similar critiques of the multivariate diagonal GARCH model have led to a plethora of alternative models of varying complexity, including the VEC, CCC, DCC, and BEKK GARCH; see FrancqZakoian2019. Analogues of these models for functional time series have not been considered, to the best of our knowledge. One potential reason for this is the fact that the identity map is not compact in general, separable Hilbert spaces; an issue that is circumvented in the formulation of pw-fARCH models. The non-compactness of the identity map prevents the user from specifying in a useful way an innovation process with identity covariance so that the conditional covariance of the model, and the model parameters, can be identified. This identifiability issue can be overcome, however, by specifying that the innovations possess some known, injective covariance operator, which is a critical idea underlying our model.

In this paper, we propose a new ARCH model, termed “operator-level ARCH”, in general separable Hilbert spaces. It provides a more direct generalization of ARCH processes to the functional setting in that the model is formulated directly for the conditional covariance operator, rather than the pointwise variance as with pw-fARCH models. Since the general model is, from a theoretical point of view, quite challenging, we focus on an operator-level counterpart of the Constant Conditional Correlation (CCC) ARCH model in the multivariate setting, which is more feasible. We establish conditions on the model parameters and innovations such that it admits strictly and weakly stationary solutions, refining the classical notion of the (functional) top Lyapunov exponent for strict stationarity to suit our general setting. Additionally, we provide sufficient conditions for the existence of finite moments and weak dependence. Our estimation approach targets the ARCH operators under known orders and relies on pseudo Yule--Walker (YW) equations. We derive consistent YW estimators for the intercept term $\Delta,$ which is a self-adjoint and positive definite operator in our setting, and the operators $\alpha_i,$ which map between operator spaces. The estimates are developed under a specific diagonal representation in finite- and infinite-dimensional settings. The identifiability issue that arises when directly deriving the YW-type equations are addressed using specific transformations of the data. In the infinite-dimensional case, we introduce a Sobolev-type condition inspired by hall:meister:2007 to quantify the approximability of infinite-dimensional operators by their finite-dimensional approximations. The estimation is conducted via Tikhonov-type pseudoinverses in the YW framework. While the theoretical properties of the estimators are well established, their complexity may present practical challenges. The applicability of the proposed methodology is illustrated through a simulation study an application to intraday log-returns.

The structure of the paper is as follows. Section (ref) introduces the general model and outlines the main assumptions. Section (ref) introduces the corresponding CCC model and derives stationarity conditions and probabilistic properties. Section (ref) addresses parameter estimation. Section (ref) presents a simulation study and Section (ref) provides applications. Section (ref) concludes. Additionally, the appendix contains preliminaries (Sections (ref) and (ref)), and proofs of key results (Section (ref)).

We adopt the following notation. The identity map is denoted by $\mathbb{I}$. On a Cartesian product space $V^n$, $n \in \mathbb{N}$, we define inner product and norm by $\langle x, y \rangle = \sum_{i=1}^n \langle x_i, y_i \rangle$ and $\|x\|^2= \sum_{i=1}^n \|x_i\|^2$ for $x = (x_1, \dots, x_n)^{\!\top},~y = (y_1, \dots, y_n)^{\!\top} \in V^n$, assuming $V$ has inner product $\langle \cdot, \cdot \rangle$ and norm $\|\cdot\|$. The spaces of bounded linear, Hilbert-Schmidt (H-S), and nuclear/trace-class operators from $\mathcal{H}$ to $\mathcal{H}_\star$ are respectively denoted by $\mathcal{L}_{\mathcal{H},\mathcal{H}_\star}$, $\mathcal{S}_{\mathcal{H},\mathcal{H}_\star}$, and $\mathcal{N}_{\mathcal{H},\mathcal{H}_\star}$, with norms $\|\cdot\|_{\mathcal{L}}, \|\cdot\|_{\mathcal{S}}, \|\cdot\|_{\mathcal{N}}$ and H-S inner product $\langle \cdot, \cdot\rangle_{\mathcal{S}}$. When $\mathcal{H}, \mathcal{H}_\star$ are clear, we write $\mathcal{T}$ for $\mathcal{T}_{\mathcal{H},\mathcal{H}_\star}$ and $\mathcal{T_H} \coloneqq \mathcal{T}_{\mathcal{H},\mathcal{H}},$ where $\mathcal{T} \in \{\mathcal{L,S,N}\}.$ For $A \in \mathcal{L}$, $A^\ast$ denotes the adjoint. Write $\mathcal{T}_{\geq 0}$ ($\mathcal{T}_{> 0}$) for the self-adjoint non-negative (positive) definite elements of $\mathcal{T}$. For $x \in \mathcal{H}, y \in \mathcal{H}_\star$, the tensor product operator is $x \otimes y \coloneqq \langle x, \cdot\rangle y$, with $x^{\otimes 2} \coloneqq x \otimes x,$ and $x\!\otimes_{\mathcal{S}}\!y$ is used when $x, y$ are H-S operators. For $p \in [1,\infty)$, $L^p_\mathcal{H} = L^p_\mathcal{H}(\Omega, \mathfrak{A},\mathbb{P})$ is the space of $X \in \mathcal{H}$ with $\mathbb{E}\|X\|^p < \infty$. The cross-covariance operator is $$ \mathscr{C}_{\!X,Y} \coloneqq \mathbb{E}[X-\mathbb{E}X ]\otimes [Y-\mathbb{E}Y], \quad X\in L^2_{\mathcal{H}},~ Y \in L^2_{\mathcal{H}_\star}, $$ and the covariance operator is $\mathscr{C}_{\!X} = \mathscr{C}_{\!X,X}$, with expectation in the Bochner sense. For weakly stationary $\boldsymbol{X}=(X_k)$ and jointly weakly stationary $\boldsymbol{X}, \boldsymbol{Y}$, the lag-$h$-covariance and lag-$h$-cross-covariance operators are $\mathscr{C}^h_{\!\boldsymbol{X}} \coloneqq \mathscr{C}_{\!X_0,X_h}$ and $\mathscr{C}^h_{\!\boldsymbol{X},\boldsymbol{Y}} \coloneqq \mathscr{C}_{\!X_0,Y_h}$, with $\mathscr{C}_{\!\boldsymbol{X}} \coloneqq \mathscr{C}^0_{\!\boldsymbol{X}}$ and $\mathscr{C}_{\!\boldsymbol{X},\boldsymbol{Y}} \coloneqq \mathscr{C}^0_{\!\boldsymbol{X},\boldsymbol{Y}}$. Further, $\boldsymbol{X} = (X_k)$ is a weak white noise (WWN) if it is weakly stationary, centered, with $\mathbb{E}\|X_0\|^2 > 0$, and $\mathscr{C}^h_{\!\boldsymbol{X}} = 0$ for all $h \neq 0$; and a \emph{strong white noise} (SWN) is an i.i.d. WWN. For sequences $(b_n), (c_n) \subset \mathbb{R}$, $b_n \asymp c_n$ (as $n \to \infty$) denotes asymptotic equivalence up to constants. Unless stated otherwise, all limits are taken as $N \to \infty$.

General model and assumptions

Let $\mathcal{H}$ be a separable Hilbert space equipped with the inner product $\langle\cdot, \cdot \rangle$.

As in the pw-fARCH framework, we call the parameter $\Delta$ the intercept term, and the parameters $\alpha_1, \dots, \alpha_p$ the (op-)ARCH operators. Here, $\Delta$ is a deterministic, self-adjoint, positive definite H-S operator, and each $\alpha_i$ is a bounded and linear operator mapping H-S to H-S operators, preserving self-adjointness and non-negative definiteness. As a consequence, $\Sigma_k$ is self-adjoint and positive definite. If (ref) admits a strictly stationary solution with finite second moments, then

so that the process exhibits conditional heteroscedasticity. Note that the model remains well-defined if $\Delta$ is merely bounded linear rather than H-S, but since consistent estimation from finite samples requires compactness, we keep the slightly stronger H-S assumption. The formula (ref) highlights the standard identifiability issue in GARCH models: neither $\Sigma_k$ nor its defining parameters are uniquely determined by (ref), since unitary transformations of $\Sigma_k$ and $\mathscr{C}_{\boldsymbol{\varepsilon}}$ leave the conditional covariance unchanged. In multivariate GARCH models this is avoided by imposing $\mathscr{C}_{\boldsymbol{\varepsilon}} = \mathbb{I}$, identifying $\Sigma_k$ as the conditional covariance operator. In infinite-dimensional Hilbert spaces, however, the identity map $\mathbb{I}$ is not compact and this is not feasible. Instead, if $\mathscr{C}_{\boldsymbol{\varepsilon}}$ is any known injective covariance operator, then (ref) uniquely identifies $\Sigma_k$. Hence, throughout we assume:

When for example $\mathcal{H}=L^2[0,1]$, and $$ \mathscr{C}_{\boldsymbol{\varepsilon}}(f)(t) = \int_0^1 C_\varepsilon(t,s)f(s)ds, $$ for a covariance kernel $C_\varepsilon$, natural choices satisfying the above are Brownian motion errors, where $C_\varepsilon(t,s) = \sigma^2 \min\{t,s\}$, and Ornstein--Uhlenbeck errors, where $C_\varepsilon(t,s) = \sigma \exp(-\theta|t-s|).$ Although the conditional covariance operator does not coincide with $\Sigma_k$, it can be easily recovered from it using the known $\mathscr{C}_{\boldsymbol{\varepsilon}}$.

In pw-f(G)ARCH models of higher order, Markovian state-space forms have been used to establish stationarity Ceroveckietal2019,Kuehnert2020. Our op-ARCH$(p)$ process also admits a Markovian representation, though in a more intricate form:

Here $X^{\otimes2,[p]}_{\!k}$ collects the past $p$ “squared” functions, so

with constant part $\boldsymbol{\Delta} \coloneqq (\Delta, 0, \dots, 0)^{\!\top}\!$, $\boldsymbol{\Psi}: \mathcal{S}^p_{\ge0}\to \mathcal{S}^p_{\ge0}$ refers to the operator-valued matrix

and $\Upsilon_{\!k}:\mathcal{S}^p_{\ge0}\to\mathcal{S}^p_{\ge0}$ refers to the map defined by $\Upsilon_{\!k}(A) \coloneqq (A_1^{1/2}\varepsilon^{\otimes2}_{k}A_1^{1/2}, A_2, \dots, A_p)^\top$ for $A=(A_1,\dots,A_p)^\top\!.$ Although the representation (ref) is natural, establishing stationarity is challenging: the nonlinearity of $\Upsilon_k$ rules out the use of a top Lyapunov exponent Kingman1973,Liggett1985, and the geometric moment contraction condition of WuShao2004, as employed in Hoermannetal2013 for pw-fARCH models, does not appear applicable to the transformation $A\mapsto\Upsilon_k(\boldsymbol{\Delta}+\boldsymbol{\Psi}(A))$. Nevertheless, a strictly stationary solution can still be constructed algorithmically FrancqZakoian2019. A more tractable model that is amenable to theoretical analysis will be introduced below.

Several of the results below, including crucially methods to estimate the operators in (ref), are simplified considerably by assuming the following:

Assumptions (ref)--(ref) imply that $\Delta$ and the ranges of $\alpha_1,\dots,\alpha_p$ are diagonalizable with respect to the (known) eigenbasis of $\mathscr{C}_\varepsilon$. In finite-dimensional (G)ARCH models, this is typically ensured by taking $\mathscr{C}_\varepsilon = \mathbb{I}$. The choice of $\mathscr{C}_\varepsilon$ may thus be viewed as selecting the basis that diagonalizes the conditional covariance. In the results below, we explicitly state which arguments rely on this assumption. Moreover, to delve more deeply into the structure of the op-ARCH model, and provide comparisons to other multivariate GARCH models, we use the notation

and assume the following.

Let $(e_j)_{j=1}^\infty$ be the eigenbasis associated to the covariance operator $\mathscr{C}_{\boldsymbol{\varepsilon}}$. Assumption (ref), together with the fact that the space of H-S operators from $\mathcal{S}^p$ to $\mathcal{S}$ is a separable Hilbert space yields

where $E_{ijk}$ is the vector placing $e_j \otimes e_k$ at position $i$ and zeros elsewhere. By the definitions of $\Sigma_k$ in the op-ARCH$(p)$ equation and $\boldsymbol{\alpha}$, and with $a_{ij\ell}\coloneqq a_{ijj\ell \ell},$ this leads to the simplified representation

Each component $\alpha_i$ of $\boldsymbol{\alpha}$ maps $\mathcal{S}_{\geq 0}$ to $\mathcal{S}_{\geq 0}$, i.e. non-negative definite, self-adjoint H-S operators into itself. Thus, for any $A = (A_1, \dots, A_p)^\top \in \mathcal{S}^p_{\geq 0}$, $\boldsymbol{\alpha}(A)$ is self-adjoint and non-negative definite under mild conditions. To establish the latter, note that we need $$ \big\langle \boldsymbol{\alpha}(A) (x), x \big\rangle \;=\; \sum^p_{i=1}\sum^\infty_{j=1}\sum^\infty_{\ell=1}\,a_{ij\ell}\langle A_i(e_j), e_j\rangle \langle x, e_\ell\rangle^2 \;\geq\; 0, \quad x \in \mathcal{H}. $$ Since $\langle A_i(e_j), e_j\rangle \geq 0$ for all $i,j$, a sufficient condition for $\boldsymbol{\alpha}(A)$ to be non-negative definite is $a_{ij\ell} \geq 0$ for all $i,j,\ell.$

CCC-op-ARCH model and structure

As indicated above, we now introduce a more parsimonious model. The expansion (ref) is the most general form of the op-ARCH operators satisfying Assumption (ref). In analogy to multivariate GARCH, this is akin to a “VEC-ARCH” specification FrancqZakoian2019. Although such a model allows for a flexible serial dependence structure, it is often considered overparameterized, challenging to estimate, and theoretically difficult to analyze. A more tractable model is obtained by retaining only the diagonal terms in (ref). This section is devoted to the discussion of such a model, where Assumptions (ref)--(ref) hold.

This is called a “CCC” model since it assumes that $\boldsymbol{\alpha}$ only depends on the diagonal-terms in its expansion, while the components of the process $X_i$ remain “conditionally correlated” through their common dependence on $\mathscr{C}_\varepsilon$. We note that one might also consider a CCC model of the form (ref) where for a compact covariance operator $\mathscr{G}$, $\Sigma_k = H_k^{1/2} \mathscr{G} H_k^{1/2}$, and $H_k$ satisfies the recursion on the right-hand side of (ref). Under the commutativity Assumption (ref) and with $\boldsymbol{\alpha}$ following (ref), this reduces to the given CCC model.

Throughout this section, $(X_k)$ is assumed to be a CCC-op-ARCH$(p)$ model for some $p\in\mathbb{N}.$ One of the benefits of this model is that there exists a linear Markovian form associated with it. To be precise, we have (see also Example (ref))

where $X^{\otimes2,[p]}_{\!k,\mathrm{d}}$ is the diagonal part of $X^{\otimes2,[p]}_{\!k}$ in (ref), which is for each component defined by

Here $\boldsymbol{\Delta}_k \coloneqq (\psi_k(\Delta), 0, \dots, 0)^\top\!,$ and $\boldsymbol{\Psi}_{\!k}$ has the same form as $\boldsymbol{\Psi}$ in (ref) with each $\alpha_i$ replaced by $\psi_k\circ\alpha_i,$ where $\circ$ refers to composition, and

It should be noted that a host of other potential, simplified models starting from (ref) might be considered. For ease of presentation and due to the empirical performance of the CCC-op-ARCH$(p)$ model in our analyses below, we have chosen to focus on this case. We discuss this and avenues for future research in Section (ref).

Strict stationarity

To deduce conditions under which this model admits a stationary solution, we introduce a top Lyapunov exponent Kingman1973,Liggett1985. To this end, we define

where $\mathcal{S}_{\ge0,\mathrm{d}}\subset \mathcal{S}_{\ge0}$ denotes the subset of non-negative, self-adjoint H-S operators with the diagonal form $A=\sum_j a_j(e_j\otimes e_j)$, and $\mathcal{B}_{\mathcal{S}^p_{\ge0,\mathrm{d}}}$ the set of bounded linear operators on $\mathcal{S}_{\ge0,\mathrm{d}}.$ Although the functional $\tau$ does not define a norm, it satisfies several useful properties for our analysis. It is dominated by the operator norm, $\tau(B) \le \|B\|_{\mathcal{L}}$, compatible with the H-S norm so that $\|B(A)\|_{\mathcal{S}} \le \tau(B)\|A\|_{\mathcal{S}},$ and sub-multiplicative, $\tau(B_1 \circ B_2) \le \tau(B_1)\tau(B_2)$.

The condition (ref) is difficult to verify in practice, even with known, relatively simple, operators $\alpha_1,...,\alpha_p$. The following two sufficient conditions are more tractable.

Existence of moments, weak dependence, and weak stationarity

Under suitable conditions, CCC-op-ARCH processes possess finite moments and display weak serial dependence. In particular, they are $L^p$-$m$-approximable, a weak dependence notion for functional data initially introduced in HoermannKokoszka2010. A process $(Y_{k})_{k\in\mathbb{Z}}\subset\mathcal{H}$ is $L^p$-$m$-approximable if (a) it admits a Bernoulli shift representation, i.e. $Y_{k} = f(\varepsilon_{k}, \varepsilon_{k-1}, \ldots)$ for some i.i.d. sequence $(\varepsilon_{k})_{k\in\mathbb{Z}}$ taking values in a measurable space $\mathbb{S}$, and a measurable mapping $f: \mathbb{S}^{\mathbb{N}} \to \mathcal{H},$ and (b) with $Y^{(m)}_{k}\!\coloneqq f(\varepsilon_{k}, \ldots, \varepsilon_{k-m+1}, \varepsilon'_{k-m}, \varepsilon'_{k-m-1}, \ldots)$, where $(\varepsilon'_k)_{k \in \mathbb{Z}}$ is an independent sequence of copies of $\varepsilon_{0},$ and with $\xi_{Y,p}(m) = (\operatorname{\mathds E}\!\|Y_{0}\!- Y^{(m)}_{0}\|_\mathcal{H}^p)^{1/p}$,

According to (a), $L^p$-m-approximable processes are strictly stationary and ergodic. If (ref) holds with $p \ge 2$, then the process $Y_k$ satisfies for example the central limit theorem. This condition is satisfied by many standard models, including functional AR and linear processes HoermannKokoszka2010, and also pw-(G)ARCH models. The following result provides sufficient conditions for this property and for the existence of finite moments. In this result, the notation $\|\cdot\|_4$ appears, referring to the norm of Schatten class operators of order 4 (see Section (ref)).

Examples

We illustrate the above results with two classes of op-ARCH processes on $\mathcal{H}=L^2[0,1]$, where $(\varepsilon_k)$ denotes an i.i.d. sequence of standard Brownian motions on $[0,1]$.

Estimation in the CCC-op-ARCH model

We now turn to the estimation of the parameters $\alpha_1,...,\alpha_p,$ and $\Delta$ in the CCC-op-ARCH$(p)$. Henceforth, $\boldsymbol{X} = (X_k) \subset \mathcal{H}$ refers to a stationary CCC-op-ARCH$(p)$ process that possesses finite second moments from which we have observed a stretch of length $N$, $X_1,...,X_N$. We will assume throughout the remainder of this article that Assumption (ref) holds so that $\mathscr{C}_\varepsilon$ and $\Sigma_k$ commute. Under this assumption, $$ \mu_{\boldsymbol{\Sigma}} = \operatorname{\mathds E}(\Sigma_k), \quad \mathscr{C}_{\!\boldsymbol{X}} = \operatorname{\mathds E}(X_k^{\otimes2}) = \mu_{\boldsymbol{\Sigma}} \mathscr{C}_{\boldsymbol{\varepsilon}}, \quad \Sigma_k^{1/2} \mathscr{C}_{\boldsymbol{\varepsilon}} \Sigma_k^{1/2}= \Sigma_k \mathscr{C}_{\boldsymbol{\varepsilon}}, \quad k \in \mathbb{Z}, $$ and the op-ARCH equations imply that

If we calculate the covariances of the above with $X^{\otimes2}_{k-i} $, $1\leq i \leq p$, due to causality the covariance with the first term on the right-hand side of (ref) will vanish, and the covariances with the second term may be expressed in terms of the $\alpha$ operators. In order to shift the nuisance operator $\mathscr{C}_\varepsilon$ off of the terms containing the $\alpha$'s, we proceed first by applying on the right a Tikhonov-regularized inverse

where $\mathbb{I}$ denotes the identity map and $\vartheta_{\!N} > 0$ is an asymptotically vanishing regularization parameter, i.e., $\vartheta_{\!N} \to 0$. After this regularization, we also project onto a finite-dimensional space. Let $(a_j, e_j)$ be the eigenpairs of $\mathscr{C}_{\boldsymbol{\varepsilon}},$ i.e. $a_1\geq a_2 \geq \cdots >0$ are the eigenvalues associated with the eigenfunctions $e_1, e_2, \dots$ of $\mathscr{C}_{\boldsymbol{\varepsilon}}.$ Then, the projection operator onto the linear space spanned by $e_1, \dots, e_K,$ $K \in \mathbb{N}$, is denoted by $\coprod^{e_K}_{e_1},$ and we define $\mathscr{C}^\ddagger_{\!\boldsymbol{\varepsilon}} \coloneqq \mathscr{C}_{\boldsymbol{\varepsilon}} \mathscr{C}^\dagger_{\!\boldsymbol{\varepsilon}}.$ According again to (ref),

Further, let $\boldsymbol{X}^{\otimes2}_{\!\boldsymbol{\varepsilon}} =(X^{\otimes2}_{\!k,\boldsymbol{\varepsilon}})_{k\in\mathbb{Z}}\subset\mathcal{S}$ and $\boldsymbol{X}^{\otimes2,[p]}\! = (X^{\otimes2,[p]}_{\!k})_{k\in\mathbb{Z}}\subset\mathcal{S}^{p}$ be the processes defined by

and their centered versions by

We see according to (ref) that the lag-1 cross-covariance operators $\mathscr{D} = \mathscr{D}_{K,N} = \mathscr{C}^1_{\!\boldsymbol{X}^{\otimes2,[p]}, \boldsymbol{X}^{\otimes2}_{\!\boldsymbol{\varepsilon}}} \in \mathcal{N}_{\mathcal{S}^p\!,\mathcal{S}}$ satisfy the following Yule–Walker (YW)-type equation:

where

with $\boldsymbol{\alpha}$ from (ref), and where the remainder $\mathscr{R} = \mathscr{R}_{K,N}$ is defined by

Notice that all the operators in the YW-type equation (ref) are well-defined due to causality of the involved processes and Proposition (ref).

If the remainder term is small, it seems natural in view of (ref) to estimate $\boldsymbol{\alpha}$ by $\hat{\mathscr{D}}\hat{\mathscr{C}}^{-1}$. This, however, does not yield a useful estimate, as $\boldsymbol{\alpha}$ is not identifiable from $\boldsymbol{\alpha}\mathscr{C}$ when $\mathscr{C}$ is not injective. For example, in $\mathcal{H} = L^2[0,1]$ with $p=1$, the operator $J \in \mathcal{S}$ with kernel $j(s,t) = -1$ if $s \leq t$ and $j(s,t) = 1$ otherwise, satisfies $\langle \mathscr{C}(K), K \rangle_{\mathcal{S}} = 0$. This issue parallels the multivariate case where, for $X \in \mathbb{R}^d$, $\mbox{vec}(XX^\top) \in \mathbb{R}^{d^2}$ does not have a full-rank covariance matrix, while the “half-vectorization” $\mbox{vech}(XX^\top) \in \mathbb{R}^{d(d+1)/2}$ typically does. In what follows, we derive estimators for the CCC-op-ARCH$(p)$ operators, that is, to reiterate, $\boldsymbol{\alpha}$ in (ref) satisfies

based on the YW-type equation (ref), suitably modified to ensure identifiability, as well as estimators for the intercept term $\Delta$.

Finite-dimensional setting

Although our ultimate goal is to derive consistent estimators for the infinite-dimensional operators in (ref), we begin with estimating the CCC-op-ARCH operators $\alpha_1, \dots, \alpha_p,$ under the simplifying assumption that

for a finite integer $K$. Further, we assume that the observed curves $X_i$ are finite-dimensional, so that with $X_{k,i} = \langle X_k, e_i\rangle:$

The representation of $X_k$ yields with $X_{k,ij}\coloneqq X_{k,i}X_{k,j}:$ \[ X^{\otimes 2}_k = \sum^K_{i=1} \sum^K_{j=1}\,X_{k,ij}(e_i\otimes e_j), \] implying that each $X^{\otimes 2, [p]}_k,$ with $p\in\mathbb{N},$ is characterized by the block matrix

Throughout, the map $\operatorname{\mathrm diag} : \mathbb{R}^K \to \mathbb{R}^{K \times K}$ and its adjoint $\operatorname{\mathrm diag}^\ast : \mathbb{R}^{K \times K} \to \mathbb{R}^K$ construct a diagonal matrix from a vector and create a vector by extracting the diagonal of the input matrix, respectively. Note that $\operatorname{\mathrm diag} : \mathbb{R}^{pK} \to \mathbb{R}^{pK \times K}$ and $\operatorname{\mathrm diag}^\ast : \mathbb{R}^{pK \times K} \to \mathbb{R}^{pK}$ are also component-wise defined, i.e. for any vectors $x_1, \dots, x_p\in\mathbb{R}^{K}$ and matrices $A_1, \dots, A_p\in\mathbb{R}^{K\times K},$

From (ref) and (ref), it follows with $\tilde{X}_{m,ij} \coloneqq X_{m,ij} - \operatorname{\mathds E}(X_{m,ij})$ for any $m$:

with $\mathscr{C}_{\mathrm{d}} \coloneqq \operatorname{\mathrm diag}^\ast\mathscr{C}\operatorname{\mathrm diag} \in \mathbb{R}^{pK\times pK},$ and where $\boldsymbol{\alpha}_{\mathrm{d}} \in \mathbb{R}^{K\times pK}$ is the block matrix $$ \boldsymbol{\alpha}_{\mathrm{d}} \coloneqq \big(\operatorname{\mathrm diag}(a_{111}, a_{122}, \dots, a_{1KK})\; \cdots \;\operatorname{\mathrm diag}(a_{p11}, a_{p22}, \dots, a_{pKK})\big). $$ Therefore, due to the identity (ref), it holds that

where $\mathscr{D}_{\mathrm{d}}\coloneqq \operatorname{\mathrm diag}^\ast\!\mathscr{D}\operatorname{\mathrm diag} \in \mathbb{R}^{pK\times K}$ and $\mathscr{R}_{\mathrm{d}}\coloneqq \operatorname{\mathrm diag}^\ast\mathscr{R}\operatorname{\mathrm diag} \in \mathbb{R}^{pK \times K},$ with $\mathscr{D}, \mathscr{R}$ from (ref). Moreover, as all the coefficients of interest $a_{ijj},$ with $1\leq i \leq p,$ $1\leq j \leq K,$ are contained in $$ \tilde{\boldsymbol{\alpha}}_{\mathrm{d}} \coloneqq \operatorname{\mathrm diag}^\ast(\boldsymbol{\alpha}_{\mathrm{d}}), $$ we propose the estimator

provided $\hat{\mathscr{C}}_{\mathrm{d}}\coloneqq\operatorname{\mathrm diag}^\ast\hat{\mathscr{C}}\operatorname{\mathrm diag}$ is non-singular, and $\hat{\mathscr{D}}_{\mathrm{d}}\coloneqq \operatorname{\mathrm diag}^\ast\!\hat{\mathscr{D}}\operatorname{\mathrm diag}$, with

The following example shows that the covariance operator of an CCC-op-ARCH process can be injective in a finite-dimensional setting.

Infinite-dimensional setting

To extend the results of the previous section to the infinite-dimensional setting, the operation “$\operatorname{\mathrm diag}$” must be generalized. In order to do so, we consider for H-S operators $A \in \mathcal{S}$ their orthogonal series expansion in the basis $(e_i\otimes e_j).$ Further we let $\ell^2 = \ell^2(\mathbb{N})$ denote the space of square-summable real valued sequences $(a_i)^\infty_{i=1} \subset \mathbb{R}.$ The bounded operator $\operatorname{\mathrm diag} : \ell^2 \to \mathcal{S}$ and its adjoint $\operatorname{\mathrm diag}^\ast : \mathcal{S} \to \ell^2$ are defined by

Here, we assume that $\boldsymbol{\alpha}$ has the infinite-dimensional form (ref). The YW-type equation (ref) also holds here, with $\mathscr{D}_{\mathrm{d}} \coloneqq \operatorname{\mathrm diag}^\ast\mathscr{D}\operatorname{\mathrm diag} : (\ell^2)^p \to \ell^2$, $\mathscr{C}_{\mathrm{d}} \coloneqq \operatorname{\mathrm diag}^\ast\mathscr{C}\operatorname{\mathrm diag} : (\ell^2)^p \to (\ell^2)^p$, and $\mathscr{R}_{\mathrm{d}} \coloneqq \operatorname{\mathrm diag}^\ast\mathscr{R}\operatorname{\mathrm diag} : (\ell^2)^p \to \ell^2$. The operator $\boldsymbol{\alpha}_{\mathrm{d}} : (\ell^2)^p \to \ell^2$ is

where each $\operatorname{\mathrm diag}((a_{ijj})^\infty_{j=1})$ is the infinite-dimensional diagonal matrix of its coefficients. Then, for any $x = (x_1, \dots, x_p)^\top \in (\ell^2)^p$, with $x_i = (x_{ij})^\infty_{j=1}$, $$ \boldsymbol{\alpha}_{\mathrm{d}}(x) = \bigg(\sum^p_{i=1}a_{ijj}x_{ij}\bigg)^\infty_{j=1}\,, $$ which lies in $\ell^2$ whenever $\sup_{i,j}|a_{ijj}| < \infty$.

Unlike in the finite-dimensional case, we cannot directly estimate all coefficients of $\boldsymbol{\alpha}_{\mathrm{d}}$, i.e. $$ \tilde{\boldsymbol{\alpha}}_{\mathrm{d}} \coloneqq \operatorname{\mathrm diag}^\ast(\boldsymbol{\alpha}_{\mathrm{d}}). $$ We therefore adopt Tikhonov regularization and project onto a finite-dimensional subspace of dimension $K \in \mathbb{N}$, with $K = K_{\!N} \to \infty$ as $N \to \infty$. To estimate $\tilde{\boldsymbol{\alpha}}_{\mathrm{d}}$, we propose

where $\hat{\mathscr{D}}_{\mathrm{d}} \coloneqq \operatorname{\mathrm diag}^\ast\hat{\mathscr{D}}\operatorname{\mathrm diag}$ with $\hat{\mathscr{D}}$ in (ref), $\hat{\mathscr{C}}_{\mathrm{d}} \coloneqq \operatorname{\mathrm diag}^\ast\hat{\mathscr{C}}\operatorname{\mathrm diag}$ with $\hat{\mathscr{C}}$ in (ref), $$\hat{\mathscr{C}}^\dagger_{\mathrm{d}} \coloneqq (\hat{\mathscr{C}}_{\mathrm{d}} + \vartheta_{\!N}\mathbb{I})^{-1} \mbox{ with } \vartheta_{\!N} \to 0,$$ and where $(\hat{\lambda}_{j,\mathrm{d}}, \hat{c}_{j,\mathrm{d}})$ and $(\lambda_{j,\mathrm{d}}, c_{j,\mathrm{d}})$ are eigenpairs of $\hat{\mathscr{C}}_{\mathrm{d}}$ and $\mathscr{C}_{\mathrm{d}}$, respectively. To establish the consistency of this estimator, we impose:

We define $(\Lambda_{\ell, \mathrm{d}})_{\ell \in \mathbb{N}}$, the reciprocal eigengaps of $\mathscr{C}_{\mathrm{d}}$ as

where

Thus $\Lambda_{\ell, \mathrm{d}}$ is the reciprocal eigengap between the eigenspace of $c_{\ell,\mathrm{d}}$ and the next distinct one. By Assumption (ref) (a), $\ell_b - \ell \leq \xi$. By part (b), the empirical analogues are

Lemma (ref) and the definitions of $\hat{\mathscr{C}}_{\mathrm{d}},\mathscr{C}_{\mathrm{d}}$ yield $\hat{\Lambda}_{\ell,\mathrm{d}} = \mathrm{O}_{\operatorname{\mathds P}}(\Lambda_{\ell,\mathrm{d}})$ for $1\leq \ell \leq K+\xi$ KuehnertRiceAue2026.

To establish consistency in the H-S norm, we impose a regularity condition that governs the approximation of $\boldsymbol{\alpha}$ by its finite-dimensional projections. This condition is analogous to the Sobolev-type smoothness assumptions introduced in hall:meister:2007 for deconvolution problems.

We note that this assumption is stricter than $\boldsymbol{\alpha}$ is H-S, which is equivalent to square-summability of the coefficients. We may now state our main consistency result.

Estimation of the Intercept term

From $\mathscr{C}_{\!\boldsymbol{X}} = \mu_{\boldsymbol{\Sigma}}\mathscr{C}_{\boldsymbol{\varepsilon}}$ and Eq. (ref), it follows

where $m_p \coloneqq \operatorname{\mathds E}(X^{\otimes2,[p]}_0)$. Accordingly, we estimate $\Delta$ by

with $\hat{\mathscr{C}}_{\!\boldsymbol{X}} = \hat{m}_1$ and $\hat{m}_p$ defined in (ref).

We next state a consistency result for $\Delta.$ Since $\mathscr{C}_{\boldsymbol{\varepsilon}}$ and $\Sigma_k$ commute for each $k,$ it holds

for some non-negative, square-summable sequence $(d_i)$. Consistency of the estimation errors for $\Delta$ is also derived based on a Sobolev condition.

Simulation Study

In this section, we present the results of simulation experiments that aim to illustrate the CCC-op-ARCH$(p)$ process, and evaluate the estimation procedures detailed in Section (ref). In each of the examples below, we view $(X_k)=\{X_i(t), \; k \in \mathbb{Z}, \; t\in [0,1]\}$ as real-valued stochastic processes taking values in the Hilbert space $\mathcal{H} =L^2[0,1]$. All analysis was done on a personal laptop in the {\tt R} programming language; r. Code that may be used to reproduce the the numerical work below is available at \url{github.com/jrvanderdoes/fungarch/}.

Implementation Details

The proposed estimators require the user to specify $K$, the dimension reduction parameter, as well as the Tikhonov parameter $\vartheta_N$. Throughout, we chose $K$ according to a modified total-variation-explained (TVE) criterion. Namely, with $(e_j)$ again denoting the eigenfunctions of $\mathscr{C}_\varepsilon$, and $X_{i,k} = \sum_{j=1}^k \langle X_i, e_j \rangle e_j$, we then choose $$ K = \min\left\{ k \; : \; \frac{\sum_{i=1}^N \| X_i - X_{i,k} \|^2}{\sum_{i=1}^N \|X_i\|^2 } \le 1-\mbox{TVE} \right\}. $$ We set TVE to $0.9$ below, unless otherwise specified. In the subsequent application to high-frequency asset price data, a 90% TVE typically resulted in a $K$ in the range of 10--15.

In order to choose $\vartheta_N$, we employ 1-step ahead cross-validation. The data of length $N$ are split into training and testing sets of size $N_{train}$ and $N_{test}$. Below we use a respective $80\%$/$20\%$ split. Since the subsequent data analysis aims to forecast pointwise conditional quantiles, we chose $\vartheta_N$ in order to minimize an integrated check-loss function measuring how well $\hat{\Sigma}_i$ may be used to predict the quantiles of $X_i$. In particular, let the $\alpha$ level check-loss function $\rho_\alpha: \mathbb{R} \to [0,\infty)$ be denoted as $$ \rho_\alpha(u)=u \times\big(\alpha-\mathds{1}_{\{u<0\}}\big). $$ Note that if a CCC-op-ARCH model has Gaussian innovations, then the pointwise conditional $\alpha$ quantile of $X_i(t)$ is $$ \sqrt{{\Sigma}^{1/2}_{j}\mathscr{C}_{\boldsymbol{\varepsilon}}{\Sigma}^{1/2}_{j}(t,t)} \times \Phi^{-1}\left(\frac{\alpha}{C_\varepsilon^{1/2}(t,t)} \right). $$ For a given fitted CCC-op-ARCH model producing forecasts of the conditional covariance operator $\hat{\Sigma}_j$, viewed as a function of the parameter $\vartheta_N$ and computed with an expanding window, we then chose $\vartheta_N$ to minimize $$ \mbox{CV}_{Err}(\vartheta_N) = \frac{1}{N_{test}} \sum_{j \in \mbox{ test set }} \int_0^1 \rho_\alpha\left(\sqrt{ \hat{\Sigma}^{1/2}_{j}\mathscr{C}_{\boldsymbol{\varepsilon}}\hat{\Sigma}^{1/2}_{j}(t,t)} \times \Phi^{-1}\left(\frac{\alpha}{C_\varepsilon^{1/2}(t,t)} \right) - X_i(t) \right)\mathrm{d}t. $$ Optimizing for the Tikhonov parameter $\vartheta_N$ is somewhat computationally intensive for large $p$, $K$, and $N$. An alternative approach that is computationally efficient, although does not necessarily lead to a consistent estimator, is to use Moore--Penrose pseudoinverses moore1920, penrose1955. In this case we modify (ref) by replacing $\mathscr{C}^\dagger_{\!\boldsymbol{\varepsilon}}\!\coprod^{e_{K}}_{e_1}$ and $\hat{\mathscr{C}}^\dagger_{\mathrm{d}}$ in the definition of $\hat{\boldsymbol{\alpha}}_{\mathrm{d}}$ with

where $(\hat\xi_i,\hat\varphi_i)$ are eigenvalues and eigenfunctions of $\hat{\mathscr{C}}$. We also compared to this estimator in our simulation experiments.

Data Generation

We simulated CCC-op-ARCH$(p)$ data $(X_k)=\{X_k(t), \; k \in \mathbb{Z}, \;t\in [0,1]\}$ as stochastic processes taking values in the space $\mathcal{H} =L^2[0,1]$. We took the error covariance operator $\mathscr{C}_{\boldsymbol{\varepsilon}}$ to be a kernel-integral operator with kernel $C_\varepsilon$ corresponding to either an Ornstein-Uhlenbeck (OU) process

or a standard Brownian motion (BM)

After generating errors with the specified covariance structure, each functional data object was simulated so that

for $i=-b, \dots, N$, with $b=100$ denoting a burn-in period that is discarded. The op-ARCH operators were constructed as

with ${\bf a}_j=(a_{1,j}, \dots, a_{d,j})$ denoting a vector of scale parameters. We further set $\Delta=\mathscr{C}_\varepsilon\,$. In the simulations below, each functional data object is simulated on a grid of $r=50$ equally spaced points on the unit interval $[0,1]$. We verified in unreported simulations that increasing the value of $r$ had a negligible impact on the reported results, although taking a small value of $r$ $(r<10)$ did negatively impact the results on model estimation error.

Spaghetti-rainbow plots illustrating the sample paths of CCC-op-ARCH$(1)$ processes of length $N=200$ are shown in Figure (ref). With $\lambda_i$, $i\in\{1,2,..\}$ denoting the ordered eigenvalues of $\mathscr{C}_\varepsilon$, Figure (ref)(a) uses ${\bf a}_1=(0, 1.6/\lambda_2, 1.6/\lambda_3, 1,6/\lambda_4)$ and OU errors, and (ref)(b) uses ${\bf a}_1=(0, 1.1/\lambda_2, 1.1/\lambda_3,1.1/\lambda_4)$ with BM errors. These sample paths share some similarity with the real data studied in the following Section (ref), including periods of volatility/heteroscedasticity.

Consistency Results

We first present results on the consistency properties of the estimators proposed in Section (ref). For each setting, we generated samples with sample sizes $N$ independently $500$ times, with $N\in \{50, 100, 250, 500, 750\}$, using OU errors (ref). We considered the $\alpha_j$ as in (ref), with either ${\bf a}_j=(0.7, 0.7)$, which we term “low-dimensional”, and ${\bf a}_j =(1,1/2^2,...,1/20^2)$, which we term “high-dimensional”.

Figure (ref) and (ref) show plots of the relative absolute error between $\Delta$ and $\hat{\Delta}$

as well as for $\alpha$ and $\hat{\alpha}$,

for increasing values of $N$.

Figure (ref) illustrates the difference in estimation error from using either the Tikhonov or Moore--Penrose inverses in defining the estimator. The plots appear to confirm the asymptotic consistency of the proposed estimators in these settings. We observed that in the low-dimensional setting, the Moore--Penrose estimator tended to perform somewhat worse than the Tikhonov based estimator, whereas in the high-dimensional case the performance between the two methods was more comparable. We note that the cross-validation criteria for determining $\vartheta_{\!N}$ does not intend to optimize the normed estimation error for the $\alpha_j$, but rather attempts to improve quantile forecasting.

Figure (ref) shows plots of $e_{N,\Delta}$ and $e_{N,\alpha}$ in the low-dimensional setting for CCC-op-ARCH$(1)$ and CCC-op-ARCH$(5)$ data. We observed decreasing estimation error as a function of $N$ in each setting, including when fitting a CCC-op-ARCH$(5)$ models to data generated according to a CCC-op-ARCH$(1)$ model (see Figure (ref)(a)).

Application to Intra-Day Return Data

ARCH models are most commonly applied to model financial return data. We considered intra-day price data of the Exchange-Traded-fund SPY that tracks the S&P500 index. The returns were created by converting the stock prices into overnight cumulative intraday returns.

The specific data we considered were obtained over two 3 year periods, 2018--2020 and 2022--2024, at a resolution of one observation every $10$ minutes ($r=39$). The full sample had a size of $N=735$ (2018--2020) and $N=717$ (2022--2024). These data and the resulting OCIDR curves in 2018--2020 are illustrated in Figure (ref).

We took as the goal of this analysis to compare the proposed model along with several other methods to forecast conditional quantiles (Value-at-Risk) of the curves $R_j$, and to evaluate the goodness-of-fit of the CCC-op-ARCH model to the data. After fitting a CCC-op-ARCH$(p)$ model, the lower $\alpha$ quantile curve forecast of $R_j(t)$ is

where $\Phi^{-1}$ is the quantile function of the standard normal distribution. In addition to comparing to the pointwise “historical” quantile computed from all the previous observations, we also fit a pw-fARCH$(1)$ model as in (ref), using the estimation method of Ceroveckietal2019, and forecast the $\alpha$ quantile curve as $\hat{V}_{j,\alpha}(t)=\hat{\sigma}_i(t)\Phi^{-1}(\alpha).$

Plots of the curves $\hat{V}_{j,0.05}(\cdot)$ for several models along with the observed OCIDR curves for a particularly volatile period in the S&P 500 index including the COVID-19 Lockdown period in March, 2020, are shown in the top panel of Figure (ref). We observed that each model exhibited a certain degree of heteroscedasticity, although this was the most pronounced for the CCC-op-ARCH$(5)$ model. The bottom panels of Figure (ref) show forecasted 95% pointwise confidence sets for the OCIDR curves based on the CCC-op-ARCH$(5)$ model.

In order to assess the accuracy of these quantile curve forecasts, we split each data set into a training set based on the first two years ($N_{train}=482$, 2018--2020 and $N_{train}=465$, 2022--2024) and then forecast a third year of data ($N_{test}=253$, 2018--2020 and $N_{test}=252$, 2022--2024), which we called the test set. Each quantile function was forecasted one step ahead using an expanding window, and the average (integrated) violation rate of each model was computed as

Table (ref) provides the observed average violation rates for each model for the nominal levels $\alpha=0.01$ and $\alpha=0.05$. We noticed that the CCC-op-ARCH$(1)$, CCC-op-ARCH$(5)$, and pw-fARCH$(1)$ models exhibited reasonably accurate coverage probabilities for the $\alpha=0.05$ level, withstanding the highly volatile S&P500 (2018--2020) sample, especially relative to the historical quantile forecast. The CCC-op-ARCH$(5)$ and pw-fARCH$(1)$ model performed well at the $\alpha=0.01$ level, with CCC-op-ARCH$(5)$ performing the best overall.

To further measure the fidelity of these forecasts to the data and contrast the forecasts of the models, we also computed the average quantile curve

for each model. These curves are shown with $\alpha=0.05$ in Figure (ref) for each of the S&P 500 samples. We observed that for 2022--2024 sample, the CCC-op-ARCH models tended to estimate a larger contrast between the variance of the curves at the beginning and end of the day, especially when compared to the pw-fARCH$(1)$ model. In the 2018--2020 sample, the average curves $\bar{V}_\alpha(t)$ were somewhat flatter for each model, although in this case only the CCC-op-ARCH$(5)$ produced forecasts with approximately nominal coverage.

In order to assess the goodness-of-fit of the estimated CCC-op-ARCH models, we computed model residuals by applying a Moore--Penrose style pseudoinverse of $\hat{\Sigma}_i$ to $X_i$. Specifically, letting $$ \hat{\Sigma}_i^\dagger = \sum_{j=1}^K\,\langle \hat{\Sigma}_i( e_j), e_j \rangle^{-1}(e_j \otimes e_j), $$ we defined residual curves

Plots of these residuals computed from CCC-op-ARCH$(p)$ models with $p\in\{1,5\}$ are shown in Figure (ref). Visually CCC-op-ARCH$(1)$ residuals appeared to retain some volatility. To evaluate for remaining conditional heteroscedasticity in the residuals, we investigated for serial correlation in the sequence of squared residual curves $Y_i(\cdot)=\hat{\varepsilon}_i^2(\cdot)$. In particular, we computed the Spherical AutoCorrelation Function (SACF), $$ \tilde{\rho}_h=\frac{1}{N-p} \sum_{i=1+p}^{N-h}\left\langle \frac{ Y_i-\mu }{\|Y_i -\mu\|}, \frac{ Y_{i+h}-\mu }{\|Y_{i+h} -\mu\|}\right\rangle, $$ as introduced in yeh:rice:2023, which is a robust estimator of autocorrelation in sequences of curves. We additionally applied the white noise test of Kokoszka:Rice:Shang:2017 to the squared residual curves Kim:Kokoszka:Rice:2023.

Plots of $\tilde{\rho}_h$ as a function of $h$ are shown in Figure (ref) for the squared residual curves derived from the S&P 500 data (2018--2020). While the original squared OCIDR curves $R_i^2$ as well as the squared residuals from the CCC-op-ARCH$(1)$ model exhibit strong serial correlation, the squared residuals from the CCC-op-ARCH$(5)$ model appear reasonably uncorrelated. Table (ref) shows the $p$-values of white noise tests applied to the squared residual series, which also support the conclusion that the CCC-op-ARCH$(1)$ model does not entirely explain the observed conditional heteroscedasticity in the data, while the CCC-op-ARCH$(5)$ model appears to fit the data well. The results were similar for the other sample considered.

Discussion

This article introduces an ARCH model for processes taking values in general separable Hilbert spaces which we call operator-level ARCH (op-ARCH) models. A key advantage of this over previous functional conditional heteroscedasticity models is that it models the complete conditional covariance function, rather than only the pointwise variance. We establish sufficient conditions for strict stationarity. Weak stationarity and the existence of finite moments and weak dependence are also discussed. Consistent operator estimates are derived both in the finite- and infinite-dimensional setting via a Yule-Walker (YW) approach. An identifiability issue complicates the direct application of YW-type equations, even under Tikhonov regularization for ill-posedness. To address this, we propose a CCC-operator-level ARCH model, which permits consistent estimation via modified YW-type equations. In finite dimensions, parametric rates are achieved, while in infinite dimensions, explicit rates depending on eigenvalue decay and operator approximation are established. An example illustrating near-parametric rates in the infinite-dimensional case is given. After detailing several aspects of implementing the proposed methods, we present results of Monte-Carlo simulation experiments, which suggest that the proposed estimators indeed appear to be consistent. In an application to cumulative intra-day return curves, the CCC-op-ARCH$(5)$ model appeared to perform well in explaining/modeling the observed heteroscedasticity in the curves, and provides alternative forecasts of the daily volatility of the curves when compared to existing pointwise models.

The model may be extended to arbitrary separable Banach spaces, drawing on the estimation frameworks of RuizAlvarez2019 and DetteKokotAue2020. Further, it would be valuable to generalize our estimation procedure to more general H-S operators. Finally, extending the theory to operator-valued GARCH processes appears promising, and work in this direction is ongoing.

Our attention was focused with regards to estimation in the “diagonal situation” (ref). Although our preliminary investigations suggest that, as in the multivariate setting, the full “VEC” model in (ref) is challenging to work with, other potential models are possible. These might include analogs of the BEKK, CCC, and DCC multivariate GARCH models, see FrancqZakoian2019. We leave this as a broad avenue for future research.

\paragraph{Acknowledgements} Parts of the article were written while Sebastian K\"uhnert was employed at University of California, Davis.

\paragraph{Funding} Alexander Aue was partially supported by NSF DMS 2515821. Sebastian K\"uhnert was partially supported by TRR 391 Spatio-temporal Statistics for the Transition of Energy and Transport (Project number 520388526) funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation). Gregory Rice was partially supported by the Discovery Grant RGPIN 50503-11525 3100 105 from the Natural Science and Engineering Research Council of Canada.