Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
74,739 characters · 24 sections · 80 citation commands
{.26in}
\thispagestyle{empty}
Bitcoin, the world's first cryptocurrency introduced in 2008, has become one of the most important digital assets worldwide. Its decentralized, peer-to-peer network, underpinned by blockchain technology, ensures transparency, security, and independence from traditional banking systems. Like other cryptocurrencies, Bitcoin is traded 24 hours a day, 7 days a week.
Recently, Bitcoin prices have increased considerably and its returns remained highly volatile. The high return volatility and speculative nature of Bitcoin present both opportunities and challenges for investors. Bitcoin's price dynamics are influenced by a complex interplay of factors, including market sentiment, regulatory developments, technological advancements, macroeconomic trends, and geopolitical events. These factors introduce substantial uncertainty, making Bitcoin forecasting particularly challenging.
Several studies have employed various methods to forecast Bitcoin returns and volatility. katsiampa2017volatility applied traditional time series models, such as GARCH, to estimate Bitcoin's volatility, providing a foundational approach to understanding its price dynamics. bouri2017return explored the return-volatility relationship in the Bitcoin market, highlighting the asset's unique risk-return profile. gradojevic2023forecasting compared the performance of Feedforward deep artificial neural network(FF-D-ANN), support vector machine(SVM) and random forest(RF) in predicting daily and hourly Bitcoin returns from two linear models: random walk (RW) and ARMAX(1,1).
Since Bitcoin is traded continuously, its returns on a given day can be viewed as functions of time measured in minutes or seconds of UTC\footnote{Coordinated Universal Time} [see zhong2024 for Bitcoin returns, and e.g. AueAlexander2017FGAC, KokoszkaPiotr2017Ifta for applications to intraday returns on stocks, stock market indexes, including S$\&$P500 and ETF's]. Then, the daily return functions can be modeled as a functional time series and predicted, using the Functional Autoregressive process or order 1 (FAR(1)). aue2015prediction show that this approach is equivalent to predicting from the Karhunen-Loeve expansion of a functional time series and Vector Autoregressive of order 1 (VAR(1)) model of eigenscores obtained from a Functional Principal Component Analysis (FPCA) of the covariance matrix. The advantage of the latter approach is that it can easily be adapted to the return dynamics characterized by conditional heteroscedasticity. The extension of the FPCA-based model to conditionally heteroscedastic data is the main objective of this paper.
We introduce a new Functional Principal Component Analysis (FPCA)-based method for forecasting daily Bitcoin return functions out of sample. Our approach aims at improved interval forecasting of daily functions of returns and exploits the serial correlation of eigenscores and squared eigenscores. The novelty of our approach is that it takes into consideration the conditional heteroscedasticity of eigenscores and produces forecast intervals that are more accurate and narrower than those obtained when the conditional heteroscedasticity is disregarded. Since the functional dynamic models for serially correlated and conditionally heteroscedastic functional time series are either complicated, or even not available in the literature, we introduce the Karhunen-Loeve (KL) dynamic factor model by imposing simplifying assumptions on the Karhunen-Loeve series expansion. This allows us to bridge the conditionally heteroscedastic dynamic discrete time dynamic models with their functional counterparts. Moreover, the KL dynamic factor model offers a convenient framework for the analysis of conditionally heteroscedastic functional time series.
Our second contribution is a rolling FPCA that allows for intraday forecasting of Bitcoin returns as subsets of a daily function. We observe that Bitcoin and many other financial assets are traded 24/7, and the delimitation of the day between hours 0:0 and 24:00 of UTC is purely conventional. This motivates us to consider the FPCA of partially overlapping daily return functions that start at subsequent time points to forecast a future segment of an incomplete function of current intraday returns.
The advantage of the functional approach compared to the conventional forecasting of discrete data is that the function can be evaluated at any time and provide forecasts at any point of time during the day rather than at a fixed sampling frequency.
In our empirical study, we examine a sample of daily functions of 15 minute and hourly returns. We model and forecast the return functions and evaluate the performance of our forecasting methods. Our method of forecasting daily return functions outperforms the FPCA-based forecasting methods available in the literature in that it produces more accurate sign forecasts and narrower forecast intervals. Our proposed algorithm for intraday return forecasting based on a rolling FPCA is shown to outperform the machine learning methods in terms of forecast accuracy and sign at fixed points of time.
In the literature, the Functional Principal Component Analysis (FPCA) has been used for forecasting first by aguilera1997approximated, who proposed a principal component prediction model using functional principal component analysis for continuous time stochastic processes. The method is applicable to forecasting long segments of future functions so that the functional time interval can be divided into the associated "past", i.e. t=1,..., T and "future", i.e. t=T+1,...,T+K, say. The FPCA is then applied separately to the divided functional data relating to the "past" and "future" times, and the eigenscores and eigenfunctions are computed for the "past" and "future" independently. The forecast future values of the function are estimated from a linear combination of the estimated "future" eigenscores and "past" eigenfunctions plus the estimated future mean. Our intraday FPCA forecasting is in part inspired by this method and allows for forecasting in real time, with the "future" potentially reduced to a neighborhood of one point in time.
aue2015prediction proposed an approach to functional forecasting of complete future functions using the FPCA combined with a Vector Autoregressive (VAR(1)) model of eigenscores. This method was applied to environmental data that follow a stationary (non-trending) process of pollution concentrations, specifically focusing on forecasting half-hourly measurements of PM10 concentrations in ambient air in Graz, Austria, from October 1, 2010, to March 31, 2011. shang2020dynamic used a similar approach, except that they used univariate autoregressive-moving average ARMA(p,q) time series models to forecast each eigenscore process independently in a study of Japanese age-specific mortality rates obtained from the Japanese Mortality Database. shang2022dynamic also applied independently the univariate ARMA(p,q) time series models to forecast the eigenscores in an empirical study of implied volatility of foreign exchange rates.
Our functional return forecasts improve upon the methods of aue2015prediction, shang2020dynamic and shang2022dynamic by adjusting them to the dynamics of financial returns. In particular, we observe that the eigenscores of returns are characterized by serial correlation at short lags, which can be efficiently modeled by the Autoregressive AR(1) process rather than an ARMA(p,q). Most importantly, we evidence serial correlation in the squared eigenscores, i.e. we show that they are conditionally heteroscedastic and characterized by the ARCH effect. This characteristic is important, and by accounting for the conditional heteroscedasticity of eigenscores, we produce improved pointwise interval forecasts of return functions.
In addition, to account for lagged cross-correlations between the eigenscores and their squares, we improve the approach of aue2015prediction by replacing the VAR(1) forecast of eigenscores, by the VAR(1)-sBEKK(1,1) model that accounts for serial correlation and conditional heteroscedasticity of eigenscores in a multivariate setup. This approach provides an additional improvement in terms of interval forecast accuracy, as the pointwise VAR-sBEKK interval forecast outperforms our AR-GARCH based method of forecasting univariate time series of eigenscores separately.
The paper is organized as follows: Section 2 describes the FPCA and KL dynamic models, introduces the new forecast method of the forecast interval and compares it with the existing FPCA forecasting models. Section 3 describes the functional data on Bitcoin and presents the results of our forecasting method for daily return functions. Section 4 introduces a new method of forecasting of returns on Bitcoin from a rolling FPCA when the function of observations is incomplete, and compares it with other Bitcoin return forecast models. Section 5 concludes the paper. Additional results are provided in the Appendices A and B..
This section reviews the Functional Principal Component Analysis (FPCA), introduces the new approach to functional times series forecasting of conditionally heteroscedastic returns, and compares it with the existing methods.
In the literature [see, e.g ramsay2005principal] functional data are usually modeled as "noise-corrupted observations on a set of trajectories" that are assumed to be realizations of a smooth random function of time $X(t)$, with unknown mean shape $\mu(t)$ and covariance $c(t,s) = cov(X(t),X(s))$ fpcatheoryliu2017functional.
Let $\{X_i(t)\}_{i \in \mathbb{Z}}$ be a stationary functional time series, where $t \in I = [a,b]$ is a continuum, and $i=1,2,...$ denote the realization. Each realization can be written as a function $X_i(t, \omega)$ in both time $t \in I$ and $\omega \in \Omega$, an element of probability space. For a fixed $\omega$, and for each $i$, $X_i(t)$ is in the Hilbert space $H = l^2(I)$, equipped with the inner product $\langle x,y\rangle_{l^2} = \int_I x(t) y(t)dt$ of square-integrable real functions on a support $I \subset \mathbb{R}$ and the norm $||x||^2_{l^2} = \int_I x^2(t) dt$ loeve1978, shang2022dynamic. At a given $t$ and for each $i$, $X_i(t, \omega)$ is a stochastic variable defined on a common probability space with finite second-order moment, i.e. any $X_i(t)$ is such that $E X_i(t)^2<\infty$ $\forall i,t$, where $E(X_i(t)^2) = \int_{\Omega} X_i^2(t, \omega) dP(\omega)$. The space $L^2(\Omega, \cal{A}, P)$ of variables with finite second-order moments is another Hilbert space equipped with the inner product $<X,Y>_{L^2} = \int_{\Omega} X(\omega) Y(\omega) dP(\omega)$.
Karhunen-Lo\`eve Expansion
The Karhunen-Lo\`eve expansion can be written for each stochastic function $X_i = \{ X_i(t),$ $t \in [a,b] \}$ that is zero-mean, or demeaned, i.e. replaced by $X_i(t) - \mu(t)$, and where $\mu(t) = E X_i(t)$ is independent of index $i$ by the stationarity assumption.
The covariance function $c(t,s)$ is:
where $\mu(t)$ is the mean function. The covariance function $c(t,s)$ allows the covariance covariance operator of $X$ denoted by $K$ to be defined as:
This operator is a function from $l^2(I)$ onto $l^2(I)$ that associates the element $K(\xi) \in l^2(I)$ with any element $\xi$ of $l^2(I)$. This operator is continuous, symmetric, positive, i.e. $<\xi, K(\xi)>_{l^2(I)} \geq 0$, $\forall \xi \in l^2(I)$.
It follows from Mercer's and Karhunen-Lo\`eve Theorems that there is an orthonormal sequence $(\xi_j)$ of functions in $l^2(I)$ and a non-increasing sequence of positive numbers $(\lambda_j)$ such that the operator $K$ admits a spectral decomposition:
where $\lambda_j, j=1,2,...$ are the eigenvalues in strictly decreasing order and $\xi_j(s), j=1,2,.. $ are the normalized eigenfunctions, so that $K(\xi_j) = \lambda_j \xi_j$ with $< \xi_j, \xi_k>_{l^2} \equiv \delta_{j,k} = 1$, if $j=k$, and $0$, otherwise. Moreover, $\xi_j$ are continuous functions if the associated $\lambda_j$ are strictly positive and the equality corresponds to both a convergence in $L^2$ for the Karhunen-Loève Theorem and a uniform convergence in Mercer's Theorem.
Then, each stochastic function $X_i(t, \omega)$, $i=1,2,...$ can be expressed as a linear combination of these basis functions. Through Karhunen-Loève expansion, a stochastic function $X(t, \omega)$ can be expressed as a linear combination of orthogonal basis functions $\xi_j$ determined by the covariance function of the process.
$$ X(t,\omega ) = \mu(t) + \sum_{j=1}^\infty \beta_{j} (\omega) \xi_j(t),$$
or, by omitting $\omega$ for ease of exposition, as:
where $\beta_{j}$ are eigenscores obtained from a projection in $l^2(I)$ of $[X(t)-\mu(t)]$ in the direction of the $j$th eigenfunction $\xi_j$, i.e. $\beta_{j} = <[X(t) - \mu(t)], \xi_j(t)>_{l^2(I)}$. The eigenscores $\beta_j, \; j=1,2,.. J $ are pairwise uncorrelated random variables with zero mean and variance $\lambda_j$ benko2009. Each of them can be interpreted as the contribution of $\xi_j(t)$ to $X(t)-\mu(t)$. The functions $\xi_j$ are continuous real-valued deterministic functions on $I=[a,b]$ that are pairwise orthogonal in $l^2(I)$, i.e. in the time domain. If $X(t)$ is Gaussian, then the random variables $\beta_j$ are also Gaussian, uncorrelated, and then stochastically independent.
Alternatively, we can write:
$$X(t, \omega) = \mu(t) + \sum_{j=1}^\infty \sqrt{\lambda_j} \beta_{j}^*(\omega) \xi_j(t),$$
to see that $\beta^*_j(\omega)$ are orthonormal in $L^2(\Omega, \cal{A}, P)$, and $\xi_j$ are orthonormal in $l^2(I)$. It is therefore sometimes said that the Karhunen-Loève expansion is bi-orthogonal (bi-orthonormal).
We consider the stationary functional time series indexed by $i=1,...,N$ and the Karhunen-Loève (KL) expansion of $X_i(t)$ for each $i$:
where the stochastic eigenscores $\beta_{i,j}^*$ are zero-mean with unit variance, uncorrelated for different $j$, stationary in $i$, but can be correlated for different $i$.
The above expansion is an equivalent representation of a stationary functional process (stationary in $i$) and imposes no other restriction, except for some weak regularity conditions. It can be written as:
where $e_i(t) = \sum_{j=J+1}^{\infty} \sqrt{\lambda_j} \beta_{i,j}^* \xi_j(t)$. The truncation error function, denoted by $e_i(t)$, is stationary, with mean zero and finite variance. The tuning parameter $J$ represents the number of eigenfunctions that are preserved in the approximation, resulting in a dimension reduction shang2022dynamic.
The representation (4) or (6) can be used to define a structural model for functional time series by specifying the dynamics of the eigenscores, either $\beta_{i,j}^*$, or $\beta_{i,j} = \sqrt{\lambda}_j \beta_{i,j}^*$. More precisely, we make the following assumption:
Assumption 1:
i) There exists a true value $J_0$ of $J$.
ii) The eigenscore series $\beta_{i|J_0} = (\beta_{i,1},...,\beta_{i,J_0})$ and $\beta_{i|\bar{J}_0} = (\beta_{i,J_0+1},\beta_{i,J_0+2},...)$ are independent.
iii) The conditional distribution of $\beta_{i|J_0} = (\beta_{i,1},...,\beta_{i,J_0})$ given $\beta_{\underline{i-1|J_0}} = (\beta_{i-1|J_0},\beta_{i-2|J_0},...$ has a parametric specification with parameter $\theta$.
iv) The remaining $\beta_{i|\bar{J}_0} = (\beta_{i,j| j > J_0})$ are independent white noises.
Definition 1: The KL dynamic factor model satisfies the decomposition (6) and Assumption 1.
It is easily checked that the KL expansion appears as a dynamic factor model with $J_0$ dynamic factors $\beta_{i|J_0}, i=1,...,J_0$ and an independent and identically distributed (i.i.d.) functional noise $e_i(t)$. The noise is partly degenerate since by construction, we have:
$$< e_i, \xi_j>_{l^2(I)} = 0, \;\forall j=1,...,J_0,$$
and the support of its distribution is in the vector space orthogonal to the space of $\xi_j, j=1,...,J_0$ in $l^2(I)$.
{\it Example 1: KL Linear Dynamic Factor Model}
For ease of exposition we consider a Gaussian framework with linear dynamics of order 1. Then, the specification is:
where
where $\beta_{i|J_0}$ is a vector of length $J_0$ containing $\beta_{i,j}, j=1,...,J_0$, the autoregressive matrix $\Phi$ of dimension $J_0 \times J_0$ has eigenvalues inside the unit circle, $(\epsilon_i), i=1,2,..$ is an independent and identically distributed (i.i.d.) white noise of dimension $J_0$, $\epsilon_i \sim N(0, \Sigma), \forall i$, assumed to be independent of the functional noise $\{e_i(t), t \in [a,b]\}, i=1,...,N$. The functional noise has a functional Gaussian distribution of infinite dimension denoted by $N(0, \sigma^2 \Omega_e)$, where $\Omega_e = Id - \sum_{j=1}^{J_0} \xi_j \xi_j'$, and $Id$ denotes the Identity operator on $l_2(I)$ and $\xi_j' \Omega_e \xi_j =0, j=1,...,J_0$. The model (7)-(8) depends on the parameters of different types. Some of them are scalar or multidimensional parameters, such as $\Phi, \Sigma, \sigma^2$, $\lambda_j, j=1,...,J_0$, while $\mu, \xi_j, j=1,...,J_0$, are functional parameters.
The above model can be viewed in another way. It is equivalent to consider the functional series $X_i(t), t \in I$, or the countable set of series $\beta_{i,N} = (\beta_{i,1},...,\beta_{i,j},...)$ that are related by a "one-to-one" linear relationship. Therefore, it is equivalent to write a Gaussian autoregressive model for the $X$'s or for the $\beta$'s. In the framework of this example, we have introduced a constrained Gaussian VAR(1) specification for the $\beta$'s. Therefore, a reduced-rank, constrained functional Gaussian FAR(1) is the associated representation of the functional process.
Let us now justify this representation, as a convenient specification of the dynamics of $\beta$'s for the Functional Vector Autoregressive FAR(1) process:
where the time $t$ is omitted for ease of exposition. Under the approach of 2000far1 the errors $(u_i, i \in \mathbb{Z})$ are centered (i.i.d.) innovations in $l_2(I)$ satisfying $E ||u_i||^4 = E[\int u_i^2(t)dt]^2 < \infty$ and $\Psi$ is a bounded linear operator satisfying the condition $\int \int \psi^2(s,t) ds dt <1$ that ensures that the above equation has a strictly stationary and causal solution such that $u_i$ is independent of $X_{i-1}, X_{i-2},..$.
Suppose, $\mu=0$ and $J_0$ is given by Assumption 1 i). Then, it follows from Aue at al (2015), Appendix A.1 that for $J=J_0$ the process can be written as
where $\epsilon_{i,j} = d_{i,j} + < u_{i},\xi_j>$ with the remainder terms $d_{i,j} = \sum_{j'=J_0+1}^{\infty} < Y_{i-1}, \xi_{j'} > < \Psi( \hat{\xi}_{j'}), \hat{\xi}_j >$ being independent white noises by Assumption 1 iv).
Next, the empirical counterpart of the above inner product $\beta_i = <X_i, \hat{\xi}_j>$ is written as follows. Let \(\mathbf{\beta}_{i} = (\beta_{i,1}, \beta_{i,2}, \dots, \beta_{i,J_0})'\) be a $(J_0 \times1) $ vector of eigenscores on day $i$, and write
where $\mathbf{\Pi_1}$ is the coefficient matrix of dimension $(J_0 \times J_0)$ with eigenvalues of modulus strictly less than 1, and $\boldsymbol{\epsilon}_i$ is a multivariate $(J_0 \times 1)$ white noise process with mean zero and variance $\Sigma $. Note that the marginal variance of each eigenscore vector being equal to an Identity matrix, implies an additional constraint $\Sigma = Id - \Pi \Pi'$ on the parameters of the model.
{\it Example 2: KL Conditionally Heteroscedastic Dynamic Factor Model}
Let us consider the model of Example 1, with equation (8) replaced by
and where $H_i$ is the conditional covariance matrix of $\beta_{i|J_0}$ of dimension $J_0 \times J_0$ given the information set ${\cal F}_{i-1}$ generated by the realizations up to and including $i-1$, and $\eta_i$ is an i.i.d. vector with mean 0 and such that $E\eta \eta' = Id$.
Then, a decomposition similar to (10) can be considered under a relaxed assumption on the errors being $L^4-m$-approximable KokoszkaPiotr2017Ifta leading to $\epsilon_{i,j}$, elements of vector $\epsilon_i = (<u_i, \tilde{\xi}_{1}>,...,<u_i, \tilde{\xi}_{J_0}>)' + \mbox{white noise}$, which is a vector of linear transformations of the conditionally heteroscedastic error plus noise. Therefore it is heteroscedastic too.
It is easy to see that in this context, a complicated multivariate functional ARCH model would be needed to represent the dynamics of conditional variances and the lagged cross-effects of the squared values of error functions. To our knowledge, such a multivariate functional model is not available. In contrast, the multivariate dynamics of conditional variance can be easily specified for the eigenscores in the context of the KL factor model, as follows.
The empirical counterpart of $\beta_i = <X_i, \hat{\xi}_j>$ can be specified as the VAR(1) process in eq. (11) with weak white noise errors and, for example, the BEKK model of conditional volatility EngleRobertF.1995MSGA. The BEKK model gets computationally costly when the number of time series data increases. Therefore, we can use its variant, called the scalar BEKK (sBEKK(1,1)) DingZhuanxin2001LSCC:
where, $\mathbf{H}_i$ is the $(J_0 \times J_0)$ multivariate conditional variance-covariance matrix of $\beta_{i}$, $a,g \in R$ with $a\geq 0,g \geq 0$ are scalar parameters, and $\boldsymbol{\epsilon}_i$ is the error vector of the VAR(1) representation of eigenscores in equation ((ref)). Matrix $\mathbf{C}$ is a lower triangular matrix of parameters of dimension $J_0 \times J_0$. We assume that the parameters of the sBEKK(1,1) model satisfy the standard stationarity and non-negativity conditions. Since the marginal variance of each eigenscore vector is equal to an Identity matrix of dimension $J_0$, the parameters of this model are constraint too, as illustrated in Appendix B.
{\it Example 3: KL Conditionally Heteroscedastic Univariate Dynamic Factor Model}
The KL factor model can be simplified further, as in practice one may prefer to model the contemporaneously uncorrelated eigenscore processes using either univariate or multivariate time series models shang2020dynamic, shang2022dynamic. In addition, for conditionally heteroscadastic processes we can specify the GARCH dynamics for the associated errors.
Let us suppose that the eigenscores $\beta_{i,1},..,\beta_{i,J_0}$ are independent and follow univariate AR($p$)-GARCH($K,L$) models. Then, the dynamic of eigenscore $\beta_{i,j}$ can be written as:
where $$\epsilon_{i,j} = \sqrt{v_{i,j}} z_{i,j}, \;\; \forall i=1,..,N, \\ $$
and where $a_{l,j} $ is the autoregressive coefficient on lag $l$ of eigenscore $j$, $\epsilon_{i,j}$ is the error term of eigenscore $j$, $\varsigma_0, \varsigma_{l,j}$ and $\zeta_{k,j}$ are the coefficients of the conditional variance, and $z_{i,j}$'s are i.i.d. N(0,1) variables gourieroux2012arch. We assume that the roots of the autoregressive polynomials $1-\sum_{l = 1}^p a_{l,j} x^l=0$ are outside the unit circle $\forall j=1,..,J_0$, $\varsigma_0>0$, and the remaining coefficients of the GARCH process satisfy the standard stationarity and non-negativity conditions. Note that the marginal variance of each eigenscore being equal to 1, implies additional constraints on the parameters of the model, which are illustrated in Appendix B.
Let us assume a KL dynamic factor model with mean zero satisfying Definition 1. The theoretical out-of-sample forecast and the conditional variance are derived below in the general case that includes Examples 1,2, and 3. $$ E(X_{N+1}| \underline{X_{N}}) = E([\sum_{j=1}^{J_0} \beta_{N+1,j} \xi_j]| \underline{X_{N}}) = \sum_{j=1}^{J_0} [ E( \beta_{N+1,j}|\underline{X_{N}})] \xi_j, $$
where $\underline{X_N}$ denotes the past curves and $E( \beta_{N+1,j} \underline{ X_{N}})$ is the forecast of the random eigenscore. By definition, the prediction of $\beta_{N+1,j}$ depends on the past curves through the $J_0$ first eigenscores only. Therefore, we have:
$$E(X_{N+1}| \underline{X_{N}}) = \sum_{j=1}^{J_0} [ E( \beta_{N+1,j}|\underline{\beta_{N|J_0}})] \xi_j.$$
It is either based on $E( \beta_{N+1,j}| \underline{\beta_{N,j}})$ when each eigenscore $\beta_{N+1,j}, \; j=1,..,J_0$ is predicted from its own past based on a univariate time series model (Example 2), or $E( \beta_{N+1,j} | \underline{\beta_{N}})$ where $\beta_{N}=(\beta_{N,1},...,\beta_{N,J_0})$ when the eigenscore $j$ is predicted from its own past and the past values of the remaining $J_0-1$ eigenscores (Examples 1 and 3).
The conditional variance at horizon 1 is:
where $\omega= \sigma^2 \sum_{j=J_0+1}^{\infty} \xi_j \xi_j'$, $\sum_{j=J_0+1}^{\infty} \xi_j \xi_j' = Id - \sum_{j=1}^{J_0} \xi_j \xi_j'$ where $Id$ denotes the identity operator on $l_2(I)$ and $\xi'$ denotes the adjoint operator of $\xi$. Note that $\sum_{j=J_0+1}^{\infty} \xi_j \xi_j'$ and $\sum_{j=1}^{J_0} \xi_j \xi_j'$ are orthogonal projectors (symmetric and idempotent) in the KL dynamic factor model.
Suppose the returns on Bitcoin (with a fixed holding period) is a stationary time series of functions $X_i(t), t \in I$, each defined over one day and observed on consecutive days $i=1,...,N$. Then, the intraday Bitcoin returns can be treated as a continuous time function.
Wee approximate the expansions of stochastic functions in (4) by
The mean $\mu(t)$ can be consistently estimated from a sample of $N$ functional observations as $\hat{\mu}_N(t) = \frac{1}{N} \sum_{i=1}^N X_i(t)$. The covariance operator $K$ can be consistently estimated from $\hat{c}_N(x) = \frac{1}{N} <X_i -\hat{\mu}_N,x> (X_i - \hat{\mu}_N)$ ramsay2005principal. aue2015prediction show that under general weak dependence assumptions these estimators are $\sqrt{N}$ consistent. KokoszkaPiotr2017Ifta show the consistency and normality of the mean and covariance estimators for functional time series that are autocorrelated and conditionally heteroscedastic.
Under the static FPCA, the functional principal components are extracted from $c(t,s)$, the estimated $\hat{K}$ yielding the estimated $\hat{\lambda}_j, \hat{\xi}_j, j=1,...,J$, where $[\hat{\xi}_1(t), \hat{\xi}_2(t),...]$ are the orthogonal sample eigenfunctions obtained from $\hat{c}(s,t) = \sum_{j=1}^J \hat{\lambda}_j \hat{\xi}_j(t) \hat{\xi}_j(s)$ where $\hat{\lambda}_1 > \hat{\lambda}_2 > \cdots \geq 0$ are the sample eigenvalues of $\hat{c}(s,t)$ under the static approach. The static FPCA approach has been extended to dynamic FPCA by aue2015prediction for serially correlated functions. Under the dynamic FPCA, the functional principal components are extracted from the long-run covariance matrix defined as $C(t,s) = \sum_{l=-\infty}^{\infty} cov (X_0(s), X_l(t))$, i.e. the marginal covariance of the stationary time series of functions, to ensure the consistency of the estimates. shang2022dynamic observe that the estimation of dynamic functional principal components using long-run covariance benefits the forecast if there are temporal dependencies in the functional data. However, shang2020dynamic claim that the functional time series method with dynamic functional principal component decomposition does not always outperform that with static functional principal component decomposition.
After the estimation, the in-sample fitted functions $\hat{X}_i(t)$, $i=1,...,N$ are obtained from the sample of functional observations $X_1(t)$ $,....,X_N(t)$ as:
where $\hat{\mu}(t) = \frac{1}{N} \sum_{i=1}^N X_i(t)$ is the estimated mean function, $\hat{\beta}_{ij}$ is the $j^{th}$ estimated principal component score for the $i^{th}$ observation where $\hat{\beta}_i = < [X(t) - \hat{\mu}(t)] \hat{\xi}_i(t)>$.
In practice, the optimal choice of truncation $J$ is crucial for the estimation accuracy. The choice criteria are discussed in shang2022dynamic, p.1028 and shang2020dynamic, p.31. A commonly used approach chooses $J$ set at the minimum, allowing to reach a given fraction $\delta$ of the cumulative covariance explained by the first $J$ leading components. For example $J = argmin_{J:J\geq 1} \{\sum_{j=1}^J \hat{\lambda}_j/\sum_{j=1}^N \hat{\lambda}_j 1\hskip -3pt \mbox{l}_{\hat{\lambda}_j>0} \geq \delta \}$ , with $\delta = 0.85 $.
This variance $\sigma^2$ of $e(t)$ can be estimated from the average of estimated daily variances of $\hat{e}_i(t) = X_i(t)-\hat{X}_i(t)$ over $N$ days, where $\hat{X}_i(t)$ is the fitted value defined in eq. (17). Next, to obtain $\hat{\omega}$, we multiply $\hat{\sigma}^2$ by $diag(Id - \sum_{j=1}^{J} \hat{\xi}_j \hat{\xi}_j')$, based on the entire sample of $N$ functions to approximate it according to eq. 15. We expect the result of $Id - \sum_{j=1}^{J} \hat{\xi}_j \hat{\xi}_j'$ to be numerically very close to an identity matrix.
In practice the R-package fda can be used to perform the FPCA, as it is done in this paper. The estimation of the dynamic parameters of eigenscore models depends on the dynamic model chosen for forecasting of eigenscores, and it is discussed in the next subsection.
In general, the forecast at horizon $h$ is computed for as follows:
where $\hat{\xi}= [\hat{\xi}_1,...,\hat{\xi}_J]$ is the set of estimated functional principal components estimated from $N$ functional observations and $\hat{\beta}_{N+h|N,j} $ is the forecast of the $j^{th}$ eigenscore. The forecast variance is
The forecast of the $j^{th}$ eigenscore is obtained from a time series forecasting method and motivated by aue2015prediction and the related papers by shang2020dynamic and shang2022dynamic who documented the serial correlation in eigenscore processes $(\beta_{i,j}), i=1,2,...$.
In this case, each future eigenscore $\hat{\beta}_{(N+1)j}, \dots, \hat{\beta}_{(N+h)j}$ for $j = 1, \dots, J$ is forecast by using the observations on the past values of each $j^{th}$ eigenscore $\hat{\beta}_{1j}, \dots, \hat{\beta}_{Nj}$ independently. We estimate the AR($p$)-GARCH($K,L$) model (14) in one step by the maximum likelihood\footnote{The R package rugarch is used to estimate the AR-GARCH in this paper.} to forecast each individual eigenscore $\beta_{(N+1)j}, \; j = 1,\dots, J$ and predict each eigenscore along with its future volatility $\nu_{(N+1),j}$.
The out-of sample forecast of eigenscore $\beta_{N+1,j}$ on day $i=N+1$ and the associated forecast of conditional variance lead to the $(1-\alpha)\%$ pointwise forecast interval calculated as:
where $\hat{\nu}_{N+1,j}$ is the forecast variance of $jth$ eigenscore estimated by the AR-GARCH for the future function on day $N+1$, $\hat{\xi}^2_j(t)$ is the estimated eigenfunction and $ z_{\alpha/2}$ denotes the quantile of the standard Normal distribution.
For the joint forecast of eigenscores, we apply the sBEKK model to forecast the future conditional volatility of eigenscores (${\beta}_{(N+1)}, \; j = 1, \dots, J_0$) on day $i=N+1$ in two steps. First, we estimate the VAR model (11) and predict the vector of future scores $\hat{\beta}_{(N+1)}$. Next, we apply the sBEKK(1,1) model (13) to the residuals $\mathbf{\hat{e}}_i, \; i=1,...N$ of the VAR(1) model estimated in the first step and forecast the conditional covariance matrix of eigenscores $H_{N+1}$ one day ahead\footnote{The R package BEKKs is used to estimate the sBEKK(1,1) model in this paper.}.
The pointwise forecast intervals at level $(1-\alpha)\%$ for the returns one-day-ahead are computed from the diagonal elements of matrix $\hat{\boldsymbol{\xi}}'(t)H_{N+1}{\hat{\boldsymbol{\xi}}}(t)$ as follows
where \(\mathbf{\hat{X}}_{N+1}\) is the forecast one-day-ahead function of returns forecast from the eigenscores $\hat{\beta}_{(N+1),j}, \;j = 1,\dots, J$, estimated by VAR(1)-sBEKK(1,1) model, and \(\hat{\boldsymbol{\xi}}(t)\) contains the estimated eigenfunctions $(\hat{\xi}_1(t),\dots, \hat{\xi}_{J_0}(t))'$
The main difference between our approach compared to aue2015prediction, shang2020dynamic and shang2022dynamic is that the variance of the forecasted $j$th score is past dependent and predicted, instead of being considered constant and estimated ex-post as in step 4 of the algorithm of aue2015prediction. In addition, the parameter estimates of the dynamic model of eigenscores account for conditional heteroscedasticity. Hence, they are more efficient, yielding more accurate forecasts of eigenscores.
In practice, the discrete time series of financial returns are characterized by conditional heteroscedasticity, in addition to potential serial correlation. Then, such a time series is modeled as an autoregressive AR process with GARCH errors, i.e., the AR-GARCH process. For this reason, we expect that the FPCA of return functions on Bitcoin generates eigenscores that are serially correlated and conditionally heteroscedastic.
Regarding the method of time series forecasting of eigenscores, aue2015prediction use the Vector Autoregressive (VAR) model. It is motivated by the fact that although the scores are contemporaneously uncorrelated, there may exist cross-correlations at higher lags that are accounted for in a multivariate model. Theorem 3.1 of aue2015prediction shows that the one-step ahead forecasts of functional time series from the KL expansion with VAR-based forecasts of eigenscores are asymptotically equivalent to the predictions based on the Functional Autoregressive process of order 1 (FAR(1)) fitted to the functional time series.
shang2020dynamic and shang2022dynamic use instead independent univariate time series forecasting models for each $\beta_{N+h,j}, j=1,..., J$. Then, each future eigenscore $\hat{\beta}_{(N+1)j}, \dots, \hat{\beta}_{(N+h)j}$ for $j = 1, \dots, J$ is forecast by using the observations on the past values of each $j^{th}$ eigenscore $\hat{\beta}_{1j}, \dots, \hat{\beta}_{Nj}$ independently from the autoregressive moving average(ARMA)(p,q) model. They report that the FPCA provides better out-of-sample forecasting performance than other functional models.
However, the ARMA model does not account for the potential presence of conditional heteroskedasticity, which can impact the accuracy of parameter estimators, the forecasts and the forecast intervals. The commonly used dynamic models for conditionally heteroskedastic discrete time series data are the Autoregressive and the Generalized Autoregressive Conditional Heteroskedasticity (ARCH and GARCH) models LIU2011724. Modeling multivariate volatilities has played an important role in economics and finance studies. There are multivariate GARCH models (MGARCH) that exist in the literature, including the so-called vec-model of BollerslevTim1988ACAP. However, this model does not guarantee a positive definite conditional covariance matrix. Other types of MGARCH models, such as the constant conditional correlation (CCC) model of BollerslevTim1990MtCi and the dynamic conditional correlation (DCC) model of EngleRobert2002DCCA, have been criticized in the literature too.
In the existing litrature, the forecast interval for FPCA models is obtained from the pointwise forecast bands algorithm. It was introduced by aue2015prediction, and applied by shang2020dynamic and shang2022dynamic. It follows the steps given below, where for ease of exposition the horizon $h=1$ is considered.
Step 1: Estimate $J$ sample eigenscore vectors $(\hat{\beta_1}, \dots, \hat{\beta}_J)$ and sample eigenfunctions
$[\hat{\xi}_1(t), \dots, \hat{\xi}_J(t)]$ using all $N$ observations.
Step 2: For $k \in [J+1, \dots, N-1]$, calculate the in-sample errors of the forecast function $$\hat{e}_{k+1}(t) =X_{k+1}(t) -\hat{X}_{k+1}(t), $$ where
and $\hat{\beta}_{(k+1)j}$ are forecasted by the chosen univariate or multivariate time series model.
Step 3: For each time point $t$, define the standard deviation of the in-sample errors $\hat{e}_{k+1}(t)$ as:
$$\gamma(t) = \sqrt{\frac{\sum\hat{e}_{k+1}(t)^2}{(N-1)-(J+1)}}$$.
Step 4: Seek turning parameters $\underline{\kappa}_\alpha, \bar{\kappa}_\alpha $ such that $\alpha \times 100\%$ of the errors satisfy
$$ -\underline{\kappa}_\alpha \gamma(t) \leq\hat{\epsilon}_k(t)\leq \bar{\kappa}_\alpha \gamma(t) $$
A practical approach to determine the pointwise forecast interval for a reasonable sample of size N, is:
for all points of time $t = 1, \dots, T$. Compared with this approach, we use instead the predicted volatility of eigenscores for improved forecast intervals.
We examine return functions over one UTC day between hours 0:00 and 24:00, and computed from 15-minute and hourly data on prices from Bitsatmp, one of the largest cryptocurrency exchanges, covering the period from January 1, 2022, to December 30, 2023. The 15-minute data comprises 729 days with 96 observations per day, totaling 69984 observations of prices observed every 15 minutes and 17472 hourly prices.
Below, we examine separately the return functions computed from data sampled at 15 minutes and hourly.
Let us define the 15-minute return on day $i = 1, \dots,N$ as $ X^*_i(t) = (ln(p_t)-ln(p_{t-1})) \times 100$ for $t \in [1,..,96]$ and $N=729$. The functional return series is displayed in Figure (ref).
Figure (ref) illustrates the daily functions of 15-minute BTC returns. Following the approach of ramsay2005principal, each function is next demeaned:
\[ X_i (t) = X_i^*(t) - \frac{1}{N} \sum_{i=1}^N X_i(t), \; t=1,...,96 \]
by substracting the 15-minute average return at t=1,...,96 over $N$ days. The demeaned return functions are displayed in Figure (ref) below.
Figure (ref) shows the daily functions of demeaned 15-minute returns over $N=729$ days.
A similar approach is applied to compute the hourly returns from hourly Bitcoin prices over $N=729$ days. The hourly returns are next demeaned. Figure (ref) shows the BTC hourly return functions.
Next, each hourly return function is demeaned:
$$ X_i (t) = X_i^*(t) - \frac{1}{N} \sum_{i=1}^N X_i(t), \; t=1,...,24 $$
by subtracting the hourly average returns at $t=1,...,24$ over $N$ days. We display the demeaned process in Figure (ref) below:
We consider one-day-ahead forecasts of return functions of 15-minute and hourly returns based on $250$ previous daily functions. This approach ensures a sufficient number of in-sample errors allowing us to compare our results with the in-sample error-based forecast intervals of aue2015prediction.
We determine the optimal number $J$ of eigenfunctions by considering the minimum number of leading components that allow us to explain a given level of the cumulative proportion of variance (CPV). We follow horvath2012, who use
\[ J_{\text{CPV}} = argmin_{J: 1 \leq J \leq N} \left\{ \frac{\sum\limits_{j=1}^{J} \hat{\lambda}_j}{\sum\limits_{j=1}^{N} \hat{\lambda}_j} \geq \delta \right\}, \]
where $\hat{\lambda}_j$ represents the $J^{th}$ estimated eigenvalue and $\delta = 85 \%$\footnote{In our 15-Minute case, it takes the first 46 (J = 46) eigenvalues to explain 85 % variance. The number J equals 16 in our hourly case.} horvath2012.
Let us illustrate the autocorrelation functions (ACF) and cross ACF for the first three eigenscores in 15-minute data.
Figure (ref) shows that the autocorrelations at lag 1 are statistically significant in all $\beta_j, j = 1,\dots, 3$, and in most eigenscores of higher order (not displayed). There also exist statistically significant lagged cross-correlations between $\beta_j, \beta_k$ for $j\neq k$, which motivates the use of the VAR model for the conditional means of eigenscores.
Figure (ref) shows that the GARCH effects exist for $\beta_j, j=1,...3$ and there are also significant laggedd cross-correlations between the squares of $\beta_j, \beta_k$ for $j\neq k$ at short lags. The serial correlation in the eigenscores and their squares motivates us to use the AR-GARCH and VAR-sBEKK models of eigenscores.
The figures displaying the ACF for hourly data can be found in the Appendix.
The root mean square error (RMSE) and mean absolute error (MAE) are used to evaluate the performance of the forecast at discrete points on each future function predicted one-step ahead out of sample. The formulas are:
Root Mean Square Error (RMSE)
\[ RMSE = \sqrt{\frac{1}{n}\sum_{t=1}^n(\hat{x}_t - x_t)^2} \]
Mean Absolute Error (MAE) \[ MAE = \frac{1}{n} \sum_{t=1}^n |\hat{x}_t - x_t| \]
where $n$ is the number of discrete forecast values in a 1-day-ahead function.
To evaluate the interval forecast accuracy, we use the interval score introduced by GneitingTilmann2007SPSR and GneitingTilmann2014PF.
For pointwise forecast intervals at the $100(1-\omega)\%$ nominal coverage probability, the lower and upper bounds at $\omega/2$ and $1-\omega/2$, denoted by $\hat{X}_{N+1}^{lb}(t)$ and $\hat{X}_{N+1}^{ub}(t)$. The score for the forecast interval is :
where $\mathbf{1}\{\cdot\}$ represents the binary indicator function, and $\omega$ is the significance level GneitingTilmann2007SPSR. The mean score for the one-day-head forecast interval is : \[ \bar{S}_\omega = \frac{1}{T}\sum_{t = 1}^T S_{\omega}[\hat{X}_{N+1}^{lb}(t), \hat{X}_{N+1}^{ub}(t), X_{N+1}(t)] \]
which rewards a narrow forecast interval if, and only if, the true future observation lies within the forecast interval shang2020dynamic.
In this section, each eigenscore $\beta_j, \; j = 1, \dots, J$ is considered an independent time series that follows the AR(1)-GARCH(1, 1) model to account for the conditional heteroskedasticity. Our approach is compared with that shang2022dynamic, where the eigenscores are modeled as ARMA(p,q) processes\footnote{The orders $p$ and $q$ are automatically selected by the auto.arima package in R.} and the conditional heteroskedasticity is disregarded .
Figures (ref) and (ref) show the one-day-ahead forecasts and forecast intervals for the 15-minute and hourly BTC return functions on 2022-12-26, respectively. The one-day-ahead forecast functions, red for AR(1)-GARCH(1,1) and green for ARMA(p,q), are calculated as one-period ahead conditional mean forecasts of AR-GARCH and ARMA(p,q) models, respectively.
The AR(1)-GARCH(1,1) forecast intervals of eigenscores are calculated using formula (ref) in Section (2.2). The ARMA(p,q) forecast intervals are calculated following aue2015prediction and described in Section (2.1).
We observe that our forecast function overlaps the forecast function of shang2020dynamic and shang2022dynamic and they are both close to the true functions of returns on Bitcoin. We also observe a difference in the forecast intervals as those obtained from the AR(1)-GARCH(1,1) model are narrower than those of aue2015prediction. This concerns both return series sampled at 15-minute and hourly intervals in Figures (ref) and (ref).
Table (ref) shows the performance statistics of one-day-ahead return function forecasts over 10 future days starting from 2022-12-26 computed from the AR(1)-GARCH(1,1) and ARMA(p,q) models of eigenscores. Columns 1 to 5 report the time series model of eigenscores (Col.1), the RMSE (Col.2), MAE (Col.3), the percentage of correctly forecast returns (Col.4) and the interval score (Col.5). The accuracy of forecast return functions evaluated at discrete points from both models is close. We applied the Diebold-Mariano test to the forecast errors of return functions at discrete points \(e_t = y_t - \hat{y}_t\), and based on the ARMA(p,q) models and AR(1)-GARCH(1,1) models of eigenscores. Both hourly (p-value of 0.362) and 15-minute (p-value of 0.201) results do not reject the null hypothesis, which suggests that the models have the same forecast accuracy level.
The difference is in the accuracy of pointwise forecast intervals. The results indicate that the forecast interval estimated by the AR(1)-GARCH(1,1) model is narrower than the forecast interval of aue2015prediction based on ARMA(p,q) as in shang2022dynamic and shang2020dynamic for both 15-minute and hourly return functions. Our method also outperforms the competitors in terms of the correct sign predictions, reported in Column 4 as "Sign", and of the score of interval predictions $\bar{S}_{\omega}$ (the lower, the better), given in Column 5. We present the forecast interval corresponding to a pointwise nominal 95% level ($\omega = 0.05$ in equation ((ref))). The actual coverage rate achieved by our method, indicating the proportion of times true future returns fall within the forecast intervals, is 98.33% for the 15-minute returns and 96.35% for the hourly returns.
In this section, we take into account the cross-correlation and conditional heteroskedasticity of eigenscores by using the VAR(1)-sBEKK(1,1) model (equation (ref)). Figures (ref) and (ref) compare the forecast function and forecast interval of 15-minute and hourly returns based on the VAR(1)-sBEKK(1,1) model (11) of eigenscores with the forecast function and forecast interval obtained from the VAR(1) model following the approach of aue2015prediction.
We observe again that the proposed method based on the VAR(1)-sBEKK(1,1) model of eigenscores produces forecast intervals that are narrower than those proposed by aue2015prediction.
Table (ref) describes the performance of pointwise forecasts and interval forecasts of the same time interval as Section 3.3.2, obtained by applying VAR(1) and VAR(1)-sBEKK(1,1) \footnote{The function "bekk_fit" from the BEKKs R package became unstable for a large number of variables includes. The functional variables of 15-minute returns are smoothed using a relatively smaller number of basis functions, so the number of eigenfunctions $J$ included in the sBEKK is small enough to keep "bekk_fit" stable. } Columns 1 to 5 report the observed frequency, followed by the RMS (Col.2), MAE (Col.3), the percentage of correct forecast sign (Col.4) and the forecast interval score (Col.5). The RMSE and MAE for our method are identical to those based on the VAR(1) model of aue2015prediction, because our estimation proceeds in two steps, and the first step produces the conditional mean forecasts. The Sign in Column 4 provides the percentage of successful sign forecast, which are also identical. The advantage of our approach is in the accuracy of the pointwise forecast interval, which can be evaluated by comparing the interval score $\mathbf{S_{BEKK}}$ of our model to the score $\mathbf{S_{AUE}}$ in Column 5 based on the VAR forecast interval of aue2015prediction. The results indicate that the VAR-sBEKK model outperforms the VAR-based forecast of aue2015prediction in terms of the accuracy of the forecast intervals in both 15-minute and hourly returns. For the 15-minute returns case, our method's coverage rate, corresponding to the nominal 95% confidence interval ($\omega = 0.05$ in equation ((ref))), is 97.91%, whereas for the hourly returns case, it stands at 98.75%.
By comparing the results in Tables 1 and 2, we find that the multivariate approach based on the VAR-sBEKK model performs the best in terms of the interval accuracy score in the 15-minute return case, while the univariate AR-GARCH model performs the best in terms of the interval accuracy score in the hourly returns.
So far, we have considered predicting the entire function of returns out-of-sample one day ahead, i.e. for $i=N+1$. Let us now consider a new approach for forecasting returns at shorter intraday horizons of $k$ time units, over 1 or more hours on a given day $i=N$ when the returns are observed only up to a given hour on that day. For assets traded 24/7, the day is conventionally determined by the UTC hours 0:00 and 24:00. A daily function of returns may as well start at a different point and cover 24 consecutive hours. This idea underlies the forecasting approach based on a "rolling FPCA" and inspired by the paper of aguilera1997approximated.
Let us consider the following scenario: On day $N$, We observe $N-1$ complete functions $X_i(t)$ on the interval \(t=[1, \dots, 24], \; i = 1,\dots, N-1\). Suppose that on day $N$ the return function is incomplete: we have observed only the returns up to and including t=23, and we wish to predict the return over the last hour, i.e in the neighborhood of $t=24$. To do that, we consider an auxiliary set of complete functions $X_i(s)^+$ on the interval $s=[0, \dots, 23]$, as illustrated in Table 3 below, observed on days $i=1,..., N$.
First, we perform the FPCA on the set of $N-1$ functions to get the eigenscores \(\hat{\beta}_{ij}, \; i = 1,\dots, N-1, \;j = 1,\dots,J\) and eigenfunctions $\hat{\xi}_l^*(t)$.
Next, we perform the FPCA on the set of $N$ auxiliary functions $X_i(s)^+, s=[0,\dots,23], \; i = 1,\dots,N$, yielding eigenscores $\hat{\alpha}_{ij}, i = 1,\dots,N, \;j = 1,\dots,J$, and eigenfunctions $\hat{\xi}_j^+(s)$. In order to predict the neighborhood of $X_N(24)$, we need to predict the eigenscores $\beta_{ij}$ associated with $X_N(t =1, \dots,24)$. This is done by a regression model.
Table (ref) provides an example of $X^+(s)$ and $X(t)$ when forecasting the return function over hour $24$ of UTC time. Note that there are 23 overlapping discrete hourly returns between the functions $X^+(s)$ and $X(t)$. For a prediction at horizon $k>1$, there will be $24-k$ overlapping discrete returns. The algorithm proceeds as follows:
Step 1: Calculate $N$ demeaned functions $X_{i}^+(s)$ by subtracting the sample mean $\bar{X}^+(s) = \frac{1}{N}\sum_{i = 1}^N X_{i}^+(s), \; s = 1-k,\dots, 24-k$.
Step 2: Calculate $N-1$ demeaned functions $X_{i}(t)$ by subtracting the sample mean $\bar{X}(t) = \frac{1}{N-1}\sum_{i = 1}^N X_{i}(t), \; t = 1,..., 24$.
Step 3: Estimate the eigenfunctions $\hat{\xi}_{j}^+(s), j = 1,\dots, J$ from the covariance operator of $X_{i}^+(s), i = 1, \dots, N$ and compute eigenscores $\hat{\alpha}_{ij} = <X_{i}^+, \hat{\xi}_{j}^+(s)>, i = 1,\dots, N; \; j = 1,\dots, J$ of $X^+(s)$ by FPCA.
Step 4: Estimate the eigenfunctions $\hat{\xi}_{l}(t), l = 1, \dots, J$ from the covariance operator of $X_{i}(t), i = 1, \dots, N-1$, and compute the eigenscores $\hat{\beta}_{il} = <X_{i}, \hat{\xi}_{l}(t)>$ by FPCA.
Step 5: For each $j=1,..,J$ run a regression model of a vector of eigenscores $\hat{\beta}_{.j}= [\beta_{1,j},...,\beta_{N-1,J}]'$ of length $N-1$ on the $(N-1) \times J$ matrix $\hat{\alpha}= [\alpha_{1,l},...,\alpha_{N-1,J}]'$, $i=1,...,N-1$:
\[ \hat{\beta}_{.j} = \hat{\alpha} c + v_j, \;\; \mbox{for} \; j = 1, \dots, J , \] where $c = [c_1,..,c_J]'$ is a vector of coefficients and $v_j$ is the regression error with mean 0 and finite variance.
Step 6: Predict the eigenscores $\beta_{N, j}, \; j = 1, \dots, J$ of length $J$ on day $N$ by using the estimated regression coefficients from Step 4 and the $1 \times J$ row vector of eigenscores $\hat{\alpha}_{N}$ estimated in Step 3 as the explanatory variable:
\[ \tilde{\beta}_{N,j} = \hat{\alpha}_{N} c \;\; \mbox{for} \; j = 1, \dots, J , \]
Step 7: The forecast over the 24th hour on the $N^{th}$ function is obtained from:
$\hat{X}_{N}(t) = \sum_{j = 1}^J \tilde{\beta}_{Nj} \hat{\xi}_{j}(t) + \bar{X}(t). $
The forecast function $\hat{X}_{N}(24)$ evaluated at $t=24$, gives the value is the forecast value of the one-hour-ahead target return.
Note that the FPCA performed on $X(t)$ and $X^+(s)$ can differ because the stationarity is satisfied over days $i=1,..N$. Hence, the covariances and their principal components can vary over the subsets of time $t$, although this variation is bounded by the square integrability of the functions.
We forecast the hourly returns on Bitcoin out-of-sample at horizon $k=1$ of one hour for a continuous 200 discrete hourly returns starting from 2019-02-07, which is the same out-of-sample period as gradojevic2023forecasting's subsample 5 for hourly BTC returns. There is no optimal way to determine how to choose the best number of daily functions to initiate the algorithm. We repeat the above process for a chosen range of daily function numbers of $[90,110]$. The results indicate that applying different numbers of daily functions only changes the MSE slightly; however, the number of daily functions applied does affect the correct forecast sign rate. Therefore, in this paper, the number of daily functions applied to the algorithm is chosen by the value that results in the best correct forecast sign rate.
In Step 5 we use various estimators of coefficients $c$. Because we expect the off-diagonal elements of the $\alpha'\alpha$ matrix to be close to zero by the orthogonality conditions, and the diagonal elements for large $j$ to be close to zero as well, we consider the Lasso and Ridge estimators, to ensure that this matrix is non-singular. The Support Vector Machine (SVM) is used as it provides optimal weighting and bias correction. We also account for potential non-linear relations between the eigenscores by using the Random forest (RF), and Neural Networks (NN).
We compare our forecast performance with the out-of-sample forecast for subsample 5 of gradojevic2023forecasting, who used the Feedforward Deep Artificial Neural Network(FF-D-ANN), SVM, RF estimators of the Random Walk (RW) and ARMAX(1,1) (Autoregressive Moving Average with exogenous inputs) models to forecast one-hour-ahead BTC returns.
Our proposed algorithm gives the best performance in terms of forecast error when the SVM method is used in Step 5, which yields an RMSE of $0.0321$ while the Ridge OLS results in the best forecast sign rate of $62.5\%$. In gradojevic2023forecasting, the Random Walk model performs best, resulting in an RMSE of $0.0363$ and a correct forecast sign rate of $50.24\%$.
So far we have discussed the forecasts at horizon $k=1$. At $K=1$, the coefficients $\alpha$ are arbitrarily close to $\beta$, so that the coefficient of determination $R^2$ of the Ridge, or Lasso regression in step 5 is arbitrarily close to 1. When the forecast horizon $k$ increases, the auxiliary set of functions $X^+(s)$ gets shifted backward in time by $k$ and the overlap between $X(t)$ and $X^+(s)$ diminishes. Then, the explanatory power of coefficients $\alpha$ decreases with the horizon $k$ and the $R^2$ of the regression in step 5 diminishes quickly. We observe that the $R^2$ is above 0.5 up to about horizon $k=10$, as illustrated in Figure 12 below.
Hence we recommend using the above intraday forecasting method at very short horizons.
In this paper, we define the KL dynamic factor model for the analysis of conditionally heteroscedastic functional processes. We also introduce two new Functional Principal Component Analysis (FPCA) based methods for forecasting intraday Bitcoin return. The first approach aims at improved interval forecasting of daily returns and exploits the serial correlation of eigenscores revealed in the FPCA literature. The novelty of our approach is that it also takes into consideration the conditional heteroscedasticity of returns and produces pointwise forecast intervals that are narrower than the existing literature. The second approach is a rolling FPCA for intraday forecasting of Bitcoin returns. It allows us to consider a sequence of partially overlapping daily return functions that start at subsequent time points and forecast the returns intradaily.
The empirical results of forecasting Bitcoin daily 15-minute and hourly return functions from the proposed methods indicate that accounting for the conditional heteroscedasticity of returns by forecasting the eigenscores using AR(1)-GARCH(1,1) and VAR(1)-sBEKK(1,1) models noticeably improves the accuracy of forecast intervals. Our proposed "rolling" FPCA method for forecasting one-step-ahead hourly Bitcoin returns shows better performance in terms of both forecast errors and forecast sign accuracy compared to other methods of Bitcoin forecasting in the literature gradojevic2023forecasting.
\nocite{KokoszkaP.2015FDFM, horman2013functional}