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.
78,702 characters · 19 sections · 61 citation commands
{.26in} \thispagestyle{empty}
Mixed causal and noncausal processes have become increasingly popular for modeling nonlinear and dual patterns in time series, especially in financial contexts such as cryptocurrency rates and commodity prices. These financial variables often exhibit bubbles and other short-lasting nonlinear and explosive patterns [hencic2015noncausal, lof2017noncausality, gourieroux2017noncausal, gourieroux2021forecast]. In this context, it is crucial to recognize that conventional measures of past dependence, often referred to as "causal dependence," may not sufficiently capture the temporal dependence of these variables. This motivates interest in mixed causal and noncausal processes, as discussed by breid1991maximum, lanne2011noncausal, lanne2013noncausal, gourieroux2016filtering, hecq2016identification, gourieroux2017local and cubadda2025sequential.\\ In the univariate framework, several studies investigated the nonlinear characteristics in financial data, especially bubbles, by using mixed causal and noncausal models [hencic2015noncausal, fries2019mixed, fries2019mixed, gourieroux2021forecast; gourieroux2021convolution; hecq2021forecasting]. Noncausal modeling techniques have also been applied to examine processes in a multivariate context (lanne2013noncausal, gourieroux2017noncausal,cubadda2019detecting, gourieroux2022nonlinear, rygh2022causal, cubadda2023detecting, cubadda2024optimization and cubadda2025sequential). Most studies, if not all, however, only consider multivariate mixed causal-noncausal models for a small number of time series and for a relatively small number of nonlinear transformations in the objective function. This paper investigates both high dimensional issues and introduces the Regularized Generalized Covariance (RGCov) estimator, an extension of the Generalized Covariance (GCov) estimator introduced in gourieroux2017noncausal,gourieroux2023generalized and further examined in jasiakneyazi and cubadda2024optimization. The GCov estimator is specifically designed for strictly stationary non-Gaussian processes, accommodating both causal and noncausal serial dependence.
The Ridge type regularization is a natural approach for the GCov. Indeed, the GCov estimator is a Portmentau type test that requires the inversion of the variance-covariance matrix. This regularization makes the eigenvalues "bigger" and prevents the matrix becoming non-invertible. A Lasso regularization would put to 0 the small-valued eigenvalues on the covariance matrix and makes it non-invertible. Note that the causal-noncausal VAR parameters are not regularized in this paper and their estimations will be used to build portfolios that are for instance free from nonlinear bubble patterns. This is different to approaches in which VARs are sparsified in a high-dimensional setting (see e.g. Hecq et al. (2023) for a postdoble selection approach of causal VARs).
The GCov estimator is favored here to the maximum likelihood estimator to identify and estimate mixed models because it does not require any distributional assumption on the errors of the process, other than their non-Gaussianity. It is a consistent, asymptotically normally distributed one-step estimator that is semiparametrically efficient. It can become parametrically efficient for well-selected nonlinear transformations. The GCov estimator involves the autocovariances of nonlinear functions of the process that summarize nonlinear and noncausal serial dependence [chan2006note]. It minimizes an objective function resembling a multivariate coefficient of determination, which includes the inverse of the variance matrix of a time series, i.e. its autocovariance at lag 0. When the dimension of the multivariate process is high and/or a high number of nonlinear functions are considered, this inversion may be difficult. This is the motivation for regularizing the variance in the formula of the GCov estimator.
The rest of the paper is as follows. In Section 2 we recap the main results about the GCov estimator. Section 3 proposes our new RGCov estimator whose behavior is investigated in Section 4 with a Monte Carlo study. Section 5 estimates a mixed causal-noncausal VAR for stocks that belong to the Rennix green index. Section 6 concludes.
\setcounter{equation}{0}
This section describes the causal-noncausal Vector Autoregressive (VAR) model and the semi-parametric GCov estimator along with the associated specification test.
Let us consider a strictly stationary process $\{Y_t\}$ of dimension $n$ satisfying a semi-parametric model:
where $g$ is a known function, $\tilde{Y}_t =(Y_t, Y_{t-1},...,Y_{T-p})$, $p$ is an integer, $(u_t)$ is an i.i.d. sequence, and $\theta$ is an unknown parameter vector. We assume that the model is well specified and the true value of the parameter $\theta$ is $\theta_0$.
An example of process (2.1) is the causal-noncausal VAR($p$) model. The multivariate causal-noncausal VAR($p$) process is:
$$Y_t = \Phi_1 Y_{t-1}+ \cdots + \Phi_p Y_{t-p} + u_t,$$
where $\theta = [vec \Phi_1', ..., vec \Phi_p']'$ \footnote{For any $m \times n$ matrix $A$ whose $j$th column is $a_j, \; j= 1, . . ., n$, $vec(A)$ denotes the column vector of dimension $mn$ defined as $vec(A) = (a_1'...,a_j',...,a_n')',$ where the prime denotes transposition.} and error $u_t$ is a multivariate non-Gaussian i.i.d. process with finite fourth-order moments. The roots of the characteristic equation $det(Id - \Phi_1 z - \cdots \Phi_p z^p) = 0$ are of a modulus either strictly greater than or smaller than one. Then, there exists a unique (strictly) stationary solution $(Y_t)$ with a two-sided representation $MA(\infty)$, which satisfies model ((ref)) with:
\centerline{$g(\tilde{Y}_t, \theta) = Y_t - \Phi_1 Y_{t-1} - \cdots - \Phi_p Y_{t-p} = u_t(\theta).$}
The causal-noncausal VAR($p$) model has been studied in gourieroux2016filtering,gourieroux2017noncausal, davis2020noncausal, rygh2022causal, cubadda2024optimization, and hall2024modelling. The error $u_t=u_t(\theta_0)$ cannot be interpreted as an innovation, because it is correlated with the past $y's$. Moreover, even though the function $g$ is linear in the current and lagged values of $Y_t$, the presence of noncausal components in $Y_t$ implies nonlinear dynamics of $Y_t$ from the calendar time perspective, with $E(Y_t|\underline{Y_{t-1}})$ nonlinear in $\underline{Y_{t-1}} = (Y_t, Y_{t-1},...)$ and conditional heteroscedasticity $V(Y_t| \underline{Y_{t-1}})$.
The presence of noncausal serial dependence can be detected through the analysis of nonlinear autocovariances, i.e. autocovariances of nonlinear functions of the observed process. This approach follows from chan2006note, who show that the presence of nonlinear dependence in strictly stationary non-Gaussian time series is revealed by the autocovariances of nonlinear transforms of that time series.
Let us consider nonlinear functions transforming a multivariate process $g(y_t;\theta)$ of dimension $n$ into a multivariate process $v_t(\theta)$ of dimension $K=Jn$ with the components $a_j[g_i(y_t;\theta)], j=1,...,J, i=1,...,n$. The transformed process:
is also serially i.i.d. when $\theta=\theta_0$. The transformed process has a dimension higher than $n$ because it is augmented by nonlinear differentiable functions of $g(Y_t;\theta)$, such as squares or logarithms, for example. Specifically, if $g (Y_t;\theta)$ has no finite fourth-order moment, then it can be replaced by a transformed multivariate process $v_t$ with a finite fourth-order moment to ensure the validity of the estimation procedure.
The advantage of the semi-parametric Generalized Covariance (GCov) estimator introduced in gourieroux2017noncausal,gourieroux2023generalized is that it does not require any distributional assumptions on the true errors $u_t=u_t(\theta_0)$, other than being i.i.d. and non-Gaussian. In addition, it is easy to compute and leads to a specification test with a known limiting distribution.
The GCov estimator of $\theta$ in model ((ref)) is defined as:
where
and $\hat{\Gamma}_T(h;\theta)$ is the sample autocovariance between $v_t(\theta)$ and $v_{t-h}(\theta)$, with $h$ denoting the lag. For a sample of $T$ observations $Y_1,...,Y_T$, the sample autocovariances are obtained from the process $\{ v_t(\theta) \}, t=1,...,T$ of dimension $K=Jn$:
for $h=0,...,H$.
The GCov estimator minimizes an objective function, which is equivalent to the sum of sample multivariate coefficients of determination computed at different lags $h$ from nonlinear transforms $v_t(\theta)$. Under the regularity conditions given in Appendix A, the GCov estimator is consistent, asymptotically normally distributed, and semi-parametrically efficient. It can achieve parametric efficiency for well-selected transformations, given in gourieroux2023generalized.
The GCov objective function can be used to test the goodness of fit of the model by testing the absence of nonlinear and linear serial correlation in the residuals: $H_0:$ $\Gamma(h) = 0$ for $h=1,2....$. The test statistic $\hat{\xi}_T^a(H) = T \sum_{h=1}^H Tr \hat{R}_T^2(h; \hat{\theta}_T)$ follows asymptotically a $\chi^2(K^2 H - (dim \theta))$ distribution where $K=Jn$, when the process $\{Y_t\}$ satisfies $g(\tilde{Y}_t; \theta_0) =u_t$ and $u_t$ is serially i.i.d.. The goodness of fit test at level $\alpha$ is conducted as follows: the null hypothesis $H_0$ is rejected when $\hat{\xi}_T (H) > \chi^2 _{1-\alpha}(K^2H - dim(\theta))$ and $H_0$ is not rejected otherwise.
For preliminary data analysis, the NLSD test of the absence of nonlinear and linear serial dependence in time series introduced in [jasiakneyazi] can be used. It is inspired by the specification test described above. The NLSD test needs to be applied prior to the estimation to determine whether there is evidence of noncausal serial dependence in the data. Thus, it corresponds to the special case of a strong white noise $Y_t=u_t$, where there is no parameter $\theta$ to be estimated. The test statistic is computed directly from the data transforms $v_t = [a_1(Y_{1,t}),...,a_1(Y_{n,t})....,a_1(Y_{1,t}),...,a_J(Y_{n,t})]'$ to test the null hypothesis of the absence of nonlinear and linear dependence based on the transformed time series of dimension $Jn=K$. It is based on the statistic: $\hat{\xi}_T(H) = T \sum_{h=1}^H Tr \hat{R}_T^2(h)$, where:
$$ \hat{R}_T^2(h)=\hat{\Gamma}_T(h) \hat{\Gamma}_T(0)^{-1} \hat{\Gamma}_T(h)' \hat{\Gamma}_T(0)^{-1}. $$
where
for $h=0,...,H$. Under the serial independence of $\{Y_t\}$, this statistic follows asymptotically a $\chi^2(K^2 H)$ distribution.
The diagonal GCov estimator $\hat{\theta}_T^d$ is obtained by minimizing the objective function with the variance matrix $\hat{\Gamma}_T^d(0; \theta)$ containing only diagonal elements of $\hat{\Gamma}_T(0; \theta)$. This estimator is described by gourieroux2017noncausal [see also cubadda2011testing] and is given by:
where
with the matrix $\hat{\Gamma}_T^d(0;\theta)$ being a diagonal matrix containing only the variances of $v_{i,t}(\theta), i=1,...,K$.
This estimator is not semiparametrically efficient, as it is not optimally weighted. Its asymptotic properties are described in gourieroux2017noncausal.
One may choose a high number of nonlinear transformations to improve the efficiency of the GCov estimator and asymptotic performance of the GCov specification test. The additional nonlinear transformation increase the dimension of matrix $\hat{\Gamma}_T(0; \theta)$, and inverting the variance matrix of a large dimension can become numerically challenging. jasiakneyazi discuss how to select the basis of transformations that can provide more information on the parameters, in addition to the linear and quadratic functions of errors $u_t$. In the causal-noncausal models with extreme risks and local explosive patterns, including the bubbles, the error does not necessarily has power moments. Moreover, some of the parameters driving those extreme risks, maybe non-identifiable from the transformations $u_t, u_t^2$ only. The following system of generators can be considered:
$$\mathcal{A} = \{ a_{t,p}(u) = |u|^p \exp(-t |u|), \; p \in \mathbb{N}, t \in [0,1] \},$$
with a countable dense subsystem given by:
$$\mathcal{A}_n = \{ a_{t_{j,n},p}(u) = |u|^p \exp(-t_{j,n}|u|), \; p \in \mathbb{N}, t_{j,n} \in [0,1], \; j=1,...,n \},$$
where $(t_{1,n},...,t_{n,n})$ is dense in $[0,1]$ when $n$ tends to infinity [see, bierens1990consistent, page 1448, for a similar approach]. The decreasing exponential functions provide square integrability of the power transforms ensuring adequate weighting.
\setcounter{equation}{0}
When $n$ is high but finite, the $ K \times K$ matrix $\hat{\Gamma}_T(0; \theta)$ is of high dimension, and its inverse may be difficult to compute numerically. This is, for example, the case of a causal-noncausal VAR process with a large number $n$ of components. In general, the dimension of the parameter $\theta=[vec\Phi_1',...,vec\Phi'_p]'$ of the VAR($p$) process is $p \times n^2$. Therefore, it is necessary to use non-linear transformations $J \geq pn$ to ensure the identifiability of the model. In the particular case of a VAR process, a large lag $p$ increases the dimension too, requiring in turn a large number of nonlinear transformations $J$. Therefore, the three parameters $n,J$ and $p$ play an important role in determining the dimension of the variance matrix. There are two issues with inverting the matrix $\hat{\Gamma}_T(0,\theta)$.
(i) The spectral decomposition of this matrix leads to close to zero smallest eigenvalues. Then, the inversion problem is in the numerical inversion given these close-to-zero eigenvalues. We refer to it as an issue of weak invertibility.
(ii) Moreover, even when the eigenvalues are sufficiently different from zero and positive, the available inversion algorithms may not be numerically efficient.
The first issue can be solved by introducing a regularized objective function, and the second one by applying an algorithm that updates the inverse of the regularized variance-covariance matrix with each consecutive observation.
To solve the problem of weak invertibility of matrix $\hat{\Gamma}_T(0; \theta)$, we introduce the Regularized (RGCov) estimator.
Definition 1: i) A Regularized (RGCov) estimator of $\theta_0$ in model ((ref)) is defined as:
where the objective function:
$$ L_T(\theta, \delta_T) = \sum_{h=1}^H Tr \hat{R}_T^2(h; \theta, \delta), $$
with
depends on the regularized variance matrix:
where $I$ is the identity matrix of dimension $K$ and $\delta>0$ is the shrinkage coefficient.
ii) More generally, a RGCov can be defined as $\hat{\theta}_T =\hat{\theta}_T(\delta_T)$, where the shrinkage parameter $\delta_T, \delta_T>0$, tends to zero, or to $\delta>0$, when $T$ tends to infinity. This includes as a special case $\delta_T=\delta>0, \forall T$, considered in i).
The regularization with fixed $\delta>0$ is used as in a ridge regression to solve the weak invertibility issue, but it leads to an estimator that is not asymptotically efficient. The asymptotic efficiency is recovered by adjusting appropriately the shrinkage parameter when the number of observations $T$ increases.
Let us consider first the general case when the sequence of positive shrinkage parameters $\delta_T \rightarrow \delta \geq 0$, and then the case when $\delta_T \rightarrow 0$.
a) In the general case, the objective function is
$$L_T(\theta, \delta_T) = \sum_{h=1}^H Tr[\hat{\Gamma}_T(h, \theta) \hat{\Gamma}_T(0, \theta, \delta_T)^{-1} \hat{\Gamma}_T(h, \theta) \hat{\Gamma}_T(0, \theta, \delta_T)^{-1}].$$
Proposition 1:
Under the regularity conditions given in Appendix A:
i) $\hat{\theta}_T=\hat{\theta}_T(\delta_T)$ is consistent of $\theta_0$
ii) $\sqrt{T} (\hat{\theta}_T - \theta_0)$ $ \stackrel{d}{\rightarrow} N(0, J(\theta_0, \delta)^{-1} I(\theta_0, \delta) J(\theta_0, \delta)^{-1})$, where matrices $J(\theta_0, \delta), I(\theta_0, \delta)$ are:
$$ J (\theta_0, \delta) = 2 \sum_{h=1}^H \left\{ \frac{\partial vec \Gamma(h; \theta_0)'}{\partial \theta} [\Gamma (0; \theta_0, \delta)^{-1} \otimes \Gamma (0; \theta_0, \delta)^{-1} ] \frac{\partial vec \Gamma(h; \theta_0)}{\partial \theta'} \right\}, $$
and
Thus the speed of convergence of the estimator and its asymptotic distribution do not depend on the way $\delta_T$ tends to $\delta$. However, its asymptotic variance-covariance matrix depends on $\delta$.
Because this estimator is not optimally weighted, it is not semi-parametrically efficient and belongs in the class of covariance estimators considered in gourieroux2017noncausal [see also cubadda2011testing].
b) In the special case $\delta=0$, the expression of the asymptotic distribution is simplified.
Proposition 2:
Under the regularity conditions given in Appendix A, if $\delta_T \rightarrow 0$ when $T \rightarrow \infty$:
i) $\hat{\theta}_T = \hat{\theta}_T(\delta_T)$ is consistent of $\theta_0$.
ii) $\sqrt{T} (\hat{\theta}_T - \theta_0) \stackrel{d}{\rightarrow} N(0, J(\theta_0)^{-1})$, where $J(\theta_0)= J(\theta_0,0)/2$ and
$$ J (\theta_0, \delta)/2 = I(\theta_0, 0)/4 = \sum_{h=1}^H \left\{ \frac{\partial vec \Gamma(h; \theta_0)'}{\partial \theta} [\Gamma (0; \theta_0)^{-1} \otimes \Gamma (0; \theta_0)^{-1} ] \frac{\partial vec \Gamma(h; \theta_0)}{\partial \theta'} \right\}. $$
where $\Gamma (0; \theta_0) = \Gamma (0; \theta_0, 0)$.
Since the asymptotic properties do not depend on the way $\delta_T$ tends to 0, these properties are the same as the properties of the GCov estimator. In particular, it is asymptotically semi-parametrically efficient [gourieroux2023generalized].
The different expressions of $J (\theta_0, \delta)$ and $I (\theta_0, \delta)$ (up to factor 4) are similar to the expressions appearing in the variance of a Generalized Least Squares (GLS) estimator in a regression model with a weighting matrix $\Omega^{-1}$ when the true weights are $\Omega_0$, say, that is:
$$(X'\Omega^{-1}X)^{-1} (X' \Omega^{-1} \Omega_0 \Omega^{-1} X) (X' \Omega^{-1} X)^{-1},$$
where $X = \frac{\partial vec \Gamma'(h, \theta_0)}{\partial \theta}$, $\Omega = \Gamma(0, \theta_0, \delta) \otimes \Gamma(0, \theta_0, \delta) $, and $\Omega_0 = \Gamma(0, \theta_0) \otimes \Gamma(0, \theta_0) $.
The regularization with $\delta_T$ can be used to build test statistics for testing the model specification, or the data for the absence of linear and nonlinear dependence.
We can use the regularized estimator $\hat{\theta}_T=\hat{\theta}_T(\delta_T)$ in the formulas of test statistics for testing the fit of the model to extend the GCov test introduced in [Gourieroux, Jasiak (2023)]). Let us define the residual-based Regularized GCov (RGCov) test statistic for testing the model specification:
Let us again first consider the general case.
Proposition 3: i) If $\delta_T \rightarrow \delta \geq 0$, when $T \rightarrow \infty$, then under the null hypothesis of independence of $u_t's$, the RGCov test statistic is asymptotically distributed as a positive combination of independent chi-square variables:
$$ \hat{\xi}_T(H, \delta_T) \stackrel{d}{\rightarrow} \sum_{l=1}^L \lambda_{l} Z_{l},$$
where $L=K^2H - dim \theta$, $Z_l, \, l=1,...,L$ are independent $Z_l \sim \chi^2(1)$ and $\lambda_{l}>0, \, \forall l$.
ii) In particular, when $dim \theta = 0$, i.e. for the Regularized NLSD (RNLSD) test statistic, we have:
$$ \hat{\xi}_T(H, \delta_T) \stackrel{d}{\rightarrow} \sum_{l=1}^{K^2} \lambda_{l}^* Z_{l}^*,$$
where $Z_l^*, \, l=1,...,K^2$ are independent $\chi^2(H)$ variables and $\lambda_{l}^*, l=1,...,K^2$ are the eigenvalues of the symmetric positive definite matrix:
$$[\Gamma(0)^{-1/2} \Gamma(0, \delta) \Gamma(0)^{-1/2}] \otimes [\Gamma(0)^{-1/2} \Gamma(0, \delta) \Gamma(0)^{-1/2}].$$
The properties of the tensor (Kronecker) product imply that $\lambda_{l}^*$ are obtained by considering all the products $\mu_j \mu_k$, where $\mu_k, \, k=1,...,K$ are the eigenvalues of $\Gamma(0)^{-1/2} \Gamma(0, \delta) \Gamma(0)^{-1/2}$.
Proof: See Appendix A.
When $dim(\theta)>0$ the test statistic (3.4) is testing the specification of the semi-parametric model, and it is applied to the residuals of the model estimated by the RGCov estimator. The case $dim(\theta)=0$ corresponds to the regularized extension of the NLSD test [jasiakneyazi]. The Regularized NLSD (RNLSD) test statistic tests the null hypothesis of the absence of linear and nonlinear dependence in the data themselves:
$$H_{0} = (\Gamma (h) = 0, \; h=1,...,H),$$ In the special case $\delta=0$, the test statistic (3.4) has a chi-square asymptotic distribution.
Proposition 4: If $\delta_T \rightarrow 0$ when $T \rightarrow \infty$, then it follows from Proposition 3 and the asymptotic equivalence of the GCov and RGCov estimators that, under the null hypothesis of independence of $u_t's$, the RGCov test statistic (3.4) has asymptotically a chi-square distribution with the degree of freedom equal to $K^2H-dim(\theta)$, for $dim(\theta) \geq 0$.
When the degree of freedom is greater than 30, than under the null hypothesis of independence of $u_t's$, the following function of the test statistic (3.4) evaluated for a given $H$ and $\delta$, and of the degree of freedom denoted by $\nu = K^2H-dim(\theta)$ with $dim(\theta) \geq 0$:
$$ \zeta(\hat{\xi}_T, \nu) = \sqrt{2 \, \hat{\xi}_T} - \sqrt{2 \nu -1 },$$
is asymptotically normally distributed:
$$ \zeta(\hat{\xi}_T, \nu) \sim N(0,1).$$
Hence, the null hypothesis can be alternatively tested using the function $\zeta$ and the asymptotically valid critical values of standard normal. This approach is particularly useful for models of large dimension $n$, or when a high number of nonlinear transformations $J$ is considered.
The implementation of the RGCov approach requires the inversion of matrices of high dimension, equal either to the total number $K$ of moments when computing $\hat{\Gamma}_T (0, \theta, \delta)$, or the number of parameters, when estimating the asymptotic variance-covariance matrix of the estimators. For numerical efficiency, the algorithms based on the Sherman-Morrison formula [sherman1949adjustment,sherman1950adjustment] can be used\footnote{The approach based on the Sherman-Morrison formula is cheaper than the computation of the inverse from either a spectral decomposition, or a Cholesky decomposition. Note also that these decompositions are not unique.}.
We first recall this formula and next explain its implementation.
Lemma 1 [sherman1949adjustment, sherman1950adjustment]: If $A$ is a symmetric positive definite matrix, then
$$(A+ xx')^{-1} = A^{-1} - \frac{A^{-1} x x' A^{-1}}{1+ x'A^{-1}x} $$
Proof: We have
$(A+xx')( A^{-1} - \frac{A^{-1} x x' A^{-1}}{1+ x'A^{-1}x}) = AA^{-1} + xx' A^{-1}(1 - \frac{1}{1+x'A^{-1}x} - \frac{x' A^{-1}x}{1+x'A^{-1}x}) = I.$
The result follows. QED
The above Lemma can be used to compute recursively the matrix $C_T=[\rho_1 I + \rho_2 \sum_{t=1}^T x_t x_t']^{-1}$ from the sample of $t=1,...,T$ observations.
Corollary 1: $$C_T = C_{T-1} - \frac{\rho_2 C_{T-1} x_Tx_T' C_{T-1}}{1+ \rho_2 x_T' C_{T-1} x_T}, \; T \geq 2.$$
Proof: We can apply Lemma 1 with $A=\rho_1 I + \rho_2 \sum_{t=1}^{T-1} x_t x_t'$ and $\sqrt{\rho_2} x_t$ for $x$. The result follows by observing that $A^{-1} = C_{T-1}$. QED
The starting value $C_1$ in this recursion is given in the next corollary:
Corollary 2:
$$C_1 = (\rho_1 I + \rho_2 x_1 x_1')^{-1} = \frac{1}{\rho_1} I - \frac{\rho_2}{\rho_1^2} \frac{x_1x_1'}{ 1+ \frac{\rho_2}{\rho_1} x_1'x_1}.$$
Proof: The result follows by Lemma 1 applied with $A=\rho_1I$ and $x= \sqrt{\rho_2} x_1$. QED
The recursive approach given in the above corollaries can be used to invert the matrix:
$$\hat{\Gamma}_T(0, \theta, \delta_T) = \delta_T I + \frac{1}{T} \sum_{t=1}^T [v_t(\theta) - \frac{1}{T}\sum_{t=1}^T v_t(\theta)][v_t(\theta) - \frac{1}{T}\sum_{t=1}^T v_t(\theta)]', $$
for given values of $T, \delta_T$ and $\theta$. The recursive formulas are then applied with $\rho_1=\delta_T$, $\rho_2 = \frac{1}{T}$, and $x_t = v_t(\theta)- \frac{1}{T}\sum_{t=1}^T v_t(\theta)$.
This inversion approach can be used when the objective function is maximized by an algorithm that requires evaluating the value of the objective function at $\hat{\theta}_T^{(j)}$, where $\hat{\theta}_T^{(j)}$ is the approximation of $\theta$ at iteration $j$ of the algorithm.
The above corollaries can also be used to compute numerically the asymptotically efficient variance-covariance matrix of the RGCov estimator when the inverse of the Hessian matrix is approximated by the inverse of the outer product of scores.
\setcounter{equation}{0}
This section evaluates the performance of the RGCov, GCov, and diagonal GCov estimators through three simulation studies, focusing on scenarios where the $K\times K$ matrix $\Gamma(0)$ is of high dimension. Specifically, in the first simulation study (Section (ref)), the high dimensionality arises from the number of variables: we set $K=30$, generated with $n=15$ and $J=2$. In the second simulation study (Section (ref)), the source of high dimensionality shifts from the number of variables to the number of linear and nonlinear transformations. In this regard, we again set $K=30$, but generated here with $n=3$ and $J=10$. Finally, in the third simulation study \ldots [continue with description].
We evaluate here the performance of the three estimators mentioned above in scenarios where the high dimensionality of $\Gamma(0,\theta)$ is driven by $n$, the number of variables. It should be noted that the performance of the RGCov estimator is here examined under the two shrinkage settings: $\delta_T = \delta$ for the case $\delta_T \to \delta$ as $T$ increases, and $\delta_T = \eta/T$ when considering the special case $\delta_T \to 0$ (see Section (ref)).
The data-generating process (DGP) is a 15-dimensional mixed VAR(1) model consisting of 14 real eigenvalues inside the unit circle, $\lambda_{1}$=$\{$-0.437, -0.374, -0.360, -0.263, -0.248, -0.201, -0.162, -0.105, -0.004, 0.320, 0.277, 0.200, 0.162, 0.164$\}$, and one real eigenvalue outside the unit circle, $\lambda_{2} = 1.5$. The error term is assumed to follow a multivariate Student-$t$ distribution with degrees of freedom $\nu=4$ and a scale matrix equal to the identity matrix.
For all three estimators, we set $H=2$, $J=2$ and define the two transformations as follows:
By combining the linear function $a_1$ with the quadratic function $a_2$, the estimator is well-equipped to capture both direct dynamics (first moment) and volatility dynamics (second moment) in the underlying data. Hence, $K = nJ = 30$ implies that matrices $\Gamma(h, \theta)$, for each $h \in \{0,1,2\}$, are of size $30 \times 30$. Therefore, inverting $\Gamma(0,\theta)$ for the GCov estimator is computationally challenging due to the curse of dimensionality. To address this issue, we investigate whether regularizing this matrix (RGCov) or considering only its diagonal elements (diagonal GCov) can lead to meaningful improvements.
We assess the performance of the three estimators by analyzing: bias, estimated variance and mean squared error (MSE). In particular, to measure the bias of the coefficient in the $i$-th row and $j$-th column of the autoregressive matrix, we use the following expression:
where $i, j = 1, \dots, n$; $\hat{\Phi}^{(m)}_{ij}$ denotes the estimated coefficient in the $i$-th row and $j$-th column of the autoregressive matrix obtained from the $m$-th replication; $\overline{\Phi}_{ij}$ represents the mean value of this coefficient across all $M$ replications and $\Phi_{ij,0}$ is its true population coefficient. To measure the estimated variance of the coefficient in the $i$-th row and $j$-th column of the autoregressive matrix, we use:
The estimated MSE for each coefficient in the matrix $\Phi$ is given by:
where $\widehat{\text{Bias}}(\hat{\Phi}_{ij})$ and $\text{Var}(\hat{\Phi}_{ij})$ are defined in (ref) and (ref), respectively. Due to the high dimensionality of $\Phi$, reporting bias, variance, and MSE for each coefficient is impractical. Instead, we summarize the performance of the estimators—GCov, diagonal GCov, and RGCov—by averaging all the coefficients of the matrix obtained from (ref), (ref), and (ref). The results are based on $M = 1000$ replications with increasing sample sizes $T = (200, 500, 800)$, and are presented in Figure (ref) for the case $\delta_T = \delta$, and in Figure (ref) for the case $\delta_T = \eta/T$. Table (ref) summarizes the results for both cases.
We begin by analyzing the case $\delta_T = \delta$, focusing on the first row of Figure (ref), which illustrates the performance of RGCov and GCov (a special case of RGCov when $\delta = 0$). In this three-dimensional plot, the $x$-axis denotes the sample size $T$, the $y$-axis represents the value of $\delta$, and the $z$-axis displays the performance metric (bias, variance, and MSE) computed by averaging the corresponding matrix entries, as described above. By examining increasing values of $\delta$, cross-validation can be applied to select the optimal $\delta$ that minimizes the respective metric. The results reveal a trade-off between the shrinkage coefficient and bias: while small values of $\delta$ (including $\delta = 0$) result in high bias due to insufficient regularization, large values of $\delta$ also increase bias in absolute terms by over-regularizing the estimator. However, between these extremes, for each $T$, there exists an optimal value of $\delta$ that minimizes the bias, indicated by gray cells in Table (ref). This is not true for the variance and MSE. Indeed, Figure (ref) amd Table (ref) show that while large shrinkage coefficients introduce bias by distorting estimates, they also stabilize the estimator by reducing variance and MSE. In other words, excessive shrinkage introduces estimation bias, but it reduces variance and MSE. Finally, the diagonal RGCov results are presented in Table (ref) and illustrated in the second row of Figure (ref). Each plot is two-dimensional—with $T$ on the $x$-axis and the analyzed metric on the $y$-axis—as $\delta$ is not involved in this estimation. The results highlight the better performance of the (R)GCov. Despite its computational simplicity, the diagonal GCov estimator loses crucial information by considering only the coefficients on the main diagonal of $\Gamma(0,\theta)$, making it less accurate and efficient than GCov and RGCov.
Let us now consider the case where the shrinkage coefficient decreases with the sample size, specifically the case $\delta_T = \eta / T$; see Table (ref) and Figure (ref). In the latter, the sample size $T$ is shown on the $x$-axis, the parameter $\eta$ on the $y$-axis, and the performance metric on the $z$-axis. Although the shrinkage coefficient, $\delta_T$, is not directly represented as an axis, this visualization is preferable to a two-dimensional graph with $\delta_T$ on the $x$-axis, as it clearly shows how the metric varies with both $\eta$ and $T$. Furthermore, placing $T$ and $\eta$ on the $x$- and $y$-axes, respectively, allows us to evaluate not only their individual effects but also their combined effect through the ratio $\eta/T$, reflecting the shrinkage coefficient definition. It should be noted that for the GCov and diagonal GCov, the results are the same as the previous case since these estimators do not depend on the shrinkage coefficient.
The results reveal a positive relationship between the shrinkage coefficient and bias. In particular, a high shrinkage coefficient leads to a negative bias; for instance, $\delta_T(\eta = 800, T = 200)$. Conversely, as the shrinkage coefficient decreases, the estimation bias shifts in the opposite direction, increasing toward positive and larger values, as highlighted by the point $\delta_T(\eta = 200, T = 800)$. It should be noted that, as the sample size increases, a larger $\eta$ is required to minimize the bias, as indicated by the gray cells in Table (ref), which mark the optimal $\eta$ for each $T$. This pattern underscores the importance of selecting a smaller $\eta$ in small samples to control bias. At the same time, larger samples require a larger $\eta$ to prevent $\delta_T$ from becoming too small and losing its regularization effect. As in the case of $\delta_T = \delta$, variance and MSE respond differently to shrinkage than bias: while large shrinkage coefficients introduce bias by distorting estimates, they also stabilize the estimator by reducing variance and MSE.
Finally, Table (ref) summarizes the frequency with which the correct mixed causal and noncausal model is identified—defined as having fourteen eigenvalues inside and one eigenvalue outside the unit circle—in the two simulation settings corresponding to $\delta_T = \delta$ and $\delta_T = \eta/T$. The results indicate that identification accuracy improves with larger sample sizes. Moreover, a relatively constant $\delta$ performs well in the first scenario, while higher values of $\eta$ become preferable as $T$ increases in the second scenario.
This section analyzes the performance of GCov, diagonal GCov, and RGCov in scenarios where the curse of dimensionality of the matrix $ \Gamma(0,\theta) $ is driven by $ J $, i.e., the number of transformations in (ref). The use of nonlinear transformations in the GCov is essential for correctly capturing the noncausal component and identifying the coefficients that result in an $i.i.d.$ error term [see gourieroux2023generalized]. However, the choice of both the number and type of nonlinear transformations is not trivial and depends on the specific characteristics of the investigated data. In other words, a limited set of nonlinear transformations, or an inadequate number of transformations, may not be sufficient to capture the underlying features of the process. Such an approach reduces the need to carefully pre-select specific transformations, providing greater flexibility in modeling diverse data characteristics. Consequently, the GCov estimator becomes more adaptable across different datasets, improving its overall reliability [see cubadda2024optimization]. However, while incorporating a large number of nonlinear transformations can improve the model’s ability to capture complex relationships, it also increases the dimensionality of the matrix $\Gamma(0,\theta)$, potentially leading to computational challenges and numerical instability. To address this issue, we investigate whether relying solely on the diagonal elements of $\Gamma(0,\theta)$ or applying the RGCov estimator to regularize the high-dimensional matrix can alleviate the adverse effects of dimensionality. This approach aims to balance the modeling flexibility—gained by including a rich set of transformations—with computational feasibility and stability.
We consider a 3-dimensional mixed VAR(1) model as DGP, characterized by two real eigenvalues inside the unit circle, $\lambda_1 = \{ 0.20, 0.41 \}$, and one real eigenvalue outside the unit circle, $ \lambda_2 = 1.5 $. As in the previous simulation study, we assume a multivariate Student-$t$ distribution with degrees of freedom $\nu = 4$ and a scale matrix equal to the identity matrix. For all three estimators, we set $H=2$ and $J=10$. In particular, we consider the following linear and nonlinear transformations: linear transformation $a_1\bigl[g_i(y_t; \theta)\bigr] = u_{i,t}$, quadratic transformation $a_2\bigl[g_i(y_t; \theta)\bigr] = u_{i,t}^2$, cubic transformation $a_3\bigl[g_i(y_t; \theta)\bigr] = u_{i,t}^3$, sign transformation $a_4\bigl[g_i(y_t; \theta)\bigr] = \text{sign}({u_{i,t}})$, absolute value transformation $ a_5\bigl[g_i(y_t; \theta)\bigr]= |u_{i,t}| $, cubic absolute transformation $ a_6\bigl[g_i(y_t; \theta)\bigr] = |u_{i,t}|^3 $, logarithmic transformation $ a_7\bigl[g_i(y_t; \theta)\bigr] = \log(|u_{i,t}|) $, squared logarithmic transformation $ a_8\bigl[g_i(y_t; \theta)\bigr] = \log(|u_{i,t}|)^2 $, cubic logarithmic transformation $a_9\bigl[g_i(y_t; \theta)\bigr] = \log(|u_{i,t}|)^3 $, and square root absolute transformation $ a_{10}\bigl[g_i(y_t; \theta)\bigr] = |u_{i,t}|^{1/2} $. These transformations are designed to capture a variety of different dynamic features in the data. The linear transformation preserves the original data structure, while the quadratic and cubic transformations capture higher-order dependencies typical in time series with volatility clustering. Transformations such as the absolute value and sign transformations are used to further model volatility and extreme fluctuations. In particular, these transformations separate dynamic volatility from other effects, such as bid-ask bounce, by highlighting the magnitude and directionality of fluctuations [see gourieroux2023generalized]. The cubic absolute transformation emphasizes large values, which are crucial for capturing extreme events, while logarithmic transformations compress large values to handle diminishing returns and improve robustness to outliers. Finally, the square root absolute transformation helps reduce the impact of large values or extreme fluctuations in the data, making it useful for stabilizing variance.
Even in this case, we assess the performance of the three estimators using bias, estimated variance, and MSE, as defined in (ref), (ref), and (ref), respectively. In particular, as discussed in Section (ref), each metric is computed element-wise and then summarized by averaging over all entries of the corresponding matrix. Figure (ref), Figure (ref), report the results for the diagonal GCov and the (R)Gcov under the settings $ \delta_T=\delta$ and $ \delta_T=\eta/T$, respectively, while Table (ref) summarize the results for all cases. Table (ref) presents the percentage of correct identifications of the true DGP, characterized by one eigenvalue outside the unit circle and two inside.
The results confirm the bad performance of the GCov estimator and reveal two key differences compared to those discussed in Section (ref). First, the diagonal GCov performs well when the high dimensionality of the matrix $\Gamma(0,\theta)$ is driven by the number of nonlinear transformations $J$, particularly in terms of variance and MSE for larger sample sizes ($T = 500, 800$). Second, a small shrinkage coefficient is required for RGCov not only to minimize bias but also to jointly optimize both bias and MSE in both scenarios: $\delta_T = \delta$ and $\delta_T = \eta/T$.
This section considers the monthly stock price series of 12 green energy companies included in the Renixx index. The Renixx index tracks the global renewable energy market, covering sectors such as wind, solar, bioenergy, geothermal, hydropower, electronic mobility, and fuel cells. It comprises 30 companies, each of which derives more than 50% of its revenues from these sectors (see \url{www.iwr.de/renixx}). Companies are selected based on the highest market capitalization of freely traded stocks (float market capitalization), ensuring that no single sector comprises more than 50% of the index. Furthermore, the combined weight of the fuel cell, mobility, and utility companies is limited to no more than 20% of the index, with no individual security exceeding a weight of 10$\%$. The Renixx is a total return performance index that reflects both price and dividend changes, quoted in euros. It is re-weighted quarterly and reviewed semi-annually, with adjustments made for new initial public offerings.
The sample we consider in this paper start in January 2008 and ends in January 2024, giving us a $T=192$ observations. We detrend each series by the cubic-spline method (see Hall and Jasiak, 2024) and divide them by their standard deviation. The plots of the price series are provided in Figure (ref). We fit mixed-VAR(1) with both GCov and RGCov estimators and compare the results. We consider four transformations of residuals up to power four of them and a number of lags included in the GCov estimator $H=10$. Under these settings, the dimension of the covariance matrix is $48\times48$, which is evidence we are in the high dimensional framework. For the RGCov estimator, we consider fixed $\delta=0.1,\dots, 0.5$ and choose the one that gives eigenvalues far from the unit root. Table (ref) demonstrates the eigenvalues of a lag polynomial of mixed-VAR(1) based on the GCov and RGCov estimators.
According to Table (ref), The GCov estimator gives eigenvalues close to one. Also, the value of the GCov test is significantly higher than the ones provided by RGCov. The estimation time is also more than twice that of RGCov. Another interesting fact is the estimates of mixed-VAR(1) based on GCov and RGCov provide different identification in terms of causal-noncausal orders. The GCov estimator provides 11 eigenvalues higher than one; however, RGCov with different regularization parameter values provides only one eigenvalue higher than one. To choose the regularization parameter of RGCov, we compare the closest eigenvalues to the unity for each of them and choose the one with a higher distance, excluding $\delta=0.1$, since it gives an eigenvalue close to zero. Based on this approach, we choose $\delta=0.3$ for more comparison with GCov for the rest of this section.
The GCov estimator in the high-dimensional setting has an invertibility issue for $\hat{\Gamma}(0)$, and using RGCov solves that problem. Table (ref) provides evidence for this fact by reporting the ten lowest eigenvalues of $\hat{\Gamma}(0)$ and $\hat{\Gamma}(0,\delta)$ for $\delta=0.3$. Let us go deeper to compare the results of GCov and RGCov with $\delta=0.3$. Figure (ref) shows the estimated residuals of each of the estimators. The residuals of RGCov((ref)) show a much more steady pattern than the residuals from GCov((ref)).
The representation theorem introduced by gourieroux2017noncausal for mixed processes distinguishes between their purely causal and noncausal latent components. For the mixed VAR($1$) model with a diagonalizable autoregressive coefficient matrix\footnote{Otherwise, we use the real Jordan form of matrix $\Phi$.} we can rewrite the matrix $\mathbf{\Phi}$ as:
where $\mathbf{J}$ is a diagonal matrix with eigenvalues of $\mathbf{\Phi}$ on the diagonal, and $\mathbf{A}$ is an invertible matrix consisting of the eigenvectors of $\mathbf{\Phi}$. Furthermore, let us denote by $n_1$ the number of eigenvalues of $\mathbf{\Phi}$ with a modulus strictly less than 1 and by $n_2$ the number of eigenvalues with a modulus strictly greater than 1, where $n_2 = n - n_1$. Then, it is possible to express $\mathbf{J}$ as follows:
where $\mathbf{J}_1$ is of dimension $n_1 \times n_1$ and has on its main diagonal the eigenvalues of $\mathbf{\Phi}$ that lie inside the unit circle, while $\mathbf{J}_2$ is of dimension $n_2 \times n_2$ and has on its diagonal the eigenvalues of $\mathbf{\Phi}$ that lie outside the unit circle. Consequently:
where $\mathbf{A}_1, \mathbf{A}_2$ represent the blocks in the decomposition of matrix $A$ as : $\mathbf{A}=(\mathbf{A}_1, \mathbf{A}_2)$, and $\mathbf{A}^1, \mathbf{A}^2$ represent the blocks in the decomposition of $\mathbf{A}^{-1}$ as $\mathbf{A}^{-1} = \left(\mathbf{A}_1, \mathbf{A}_2\right)^{\prime}$.\\ Since the eigenvalues of $\mathbf{J}_{1}$ in ((ref)) are in modulus less than 1, $\mathbf{y}_{1,t}^*$ is defined as the latent purely causal component of $\mathbf{y}_t$. By similar reasoning, $\mathbf{y}^*_{2,t}$ captures the noncausal component of the investigated process. In particular, $\mathbf{y}^*_{2,t}$ is the locally explosive component of $\mathbf{y}_t$ that follows a strictly stationary noncausal (V)AR process. Figure (ref) illustrates the causal and noncausal components of VAR(1) estimated with RGCov estimator with $\delta=0,3$.
As highlighted in hall2024modelling and following engle1987co, equation ((ref)) suggests that the causal and noncausal components are represented by linear combinations that eliminate the noncausal and causal components, respectively. Specifically, in a mixed VAR($1$) process, we have $n_1 + n_2 = n$, which implies $n_1 < n$ and $n_2 < n$, hence indicating the presence of common causal and noncausal components by definition. Conversely, if the process is either purely causal or purely noncausal, then $n = n_1$ or $n = n_2$, respectively. This means that if the process is either purely causal or purely noncausal, each series in the investigated data is driven by independent causal or noncausal components, assuming that $\mathbf{A}$ is a full-rank matrix. Hence, the causal component $\mathbf{y}^*_{1,t}$ is the stationary linear combination that eliminates its local explosive characteristics and can be interpreted as bubble 'cointegration' in the sense that it removes the nonlinear patterns captured by the noncausal component [see hall2024modelling].\footnote{Note that cubadda2023detecting and cubadda2019detecting investigate the presence of co-movements in the multiplicative VAR model, necessitating the choice of models in terms of lag and lead polynomials.} In the next subsection, we propose portfolio management based on the causal and noncausal components of the VAR(1) estimated by RGCov with $\delta=0.3$.
In the framework of causal-noncausal VAR models, bubbles arise as inherent features of a strictly stationary multivariate process. The bubbles are common to the component series, estimable, and predictable because the bubble dynamics are approximated by the non-causal component of the VAR process. In addition, the linear combinations of green stocks that eliminate common bubbles follow directly from the state-space representation of the causal-noncausal VAR model and are estimable and predictable as well from the causal components of the VAR. Hence, the causal and non-causal components, including the bubble, can be interpreted as ”common features” of the green stocks in the sense of engle1993testing and compared with the bubble cointegration of cubadda2024optimization. The causal component of the VAR model is a bubble-free linear combination of green stocks. It represents stable investment portfolios, while the portfolios that "ride the bubble" can be obtained from the noncausal combinations of green stocks. The latter portfolios assume more risk but provide higher returns during the bubble episode.
Let us consider the causal and the noncausal components of the mixed VAR(1) fitted to twelve green stocks in the Renixx index. We investigate the performance of portfolios of these stocks with the allocations determined from the coefficients of the causal and the noncausal components of the process. The objective is to compare the performance of the Rennix index with the stable portfolio and the portfolio with the allocations determined from the explosive noncausal component in terms of cumulative returns. To construct the portfolios, we focus our attention on the current causal component $A^1_t$ and noncausal component $A^2_t$.
We consider an investor who invests $V_1=100\$$ in the green stocks at time $t=1$ and sells the stocks at the end of the month for the amount of $V_2=V_1+r_2$, where $r_2$ is the return over the first month. At time $t=2$, s/he invests $V_2=V_1+r_2$. Therefore, at any time $t=1,..T$, the investor invests $V_t=V_1 + \sum_{h=2}^t r_h$ and sells the stocks at the end of each month. This strategy is repeated for $t=1,2,..., T-1$. The negative coefficients in the allocation vector provided by the causal and non-causal components are interpreted as a short sell. The initial investment of $V_1$ is equal for all of the portfolios.
We can construct 12 causal and noncausal, i.e. stable portfolios from the allocations given in $A^1_t$ and $A^2_t$. We consider $a_{ij}$ indicates $i^{th}$ portfolio and $j^{th}$ asset. According to the budget constraints on $V_t$ in each period, the weights are equal to $w_{ijt}=s_{it}* a_{ij}$ where $s_{it} $ is:
$$s_{it}= \frac{V_{it}}{\sum_{j=1}^{12} |a_{ij}|*p_{j,t}},$$
where $V_{it}$ is the value of the investment at time $t$ in portfolio $i$ and $p_{jt}$ is the price of asset j at time t. To estimate the return of each portfolio, we have $$r_{i,t+1}= \sum_{j=1}^{12} w_{ijt}*p_{j{t+1}} - \sum_{j=1}^{12} w_{ijt}*p_{j{t}}.$$ $$r_{i,t+1}= s_{it} \sum_{j=1}^{12} b_{ij} [( ln(p_{j,{t+1}}) -ln(p_{j,{t}})].$$ Consequently, the cumulative sums of returns are the summation of the returns over time.
Since Renixx is an index of green stocks with positive weights, we calculate the returns on Rennix as $$r_{Renixx,t+1}= s_{Renixx,t}*[( ln(p_{Renixx,{t+1}}) -ln(p_{Renixx,{t}})],$$ where $$s_{Renixx,t}=\frac{V_{it}}{p_{Renixx,{t}}}.$$
In Figure (ref), we provide the cumulated return of the noncausal portfolio and the best causal portfolio among 11 causal portfolios and the Renixx itself. Both proposed portfolios always eliminate Renixx in terms of cumulated returns. During the bubble period, the causal portfolio performs better. However, the noncausal portfolio with the bubble effect performance is steadier overall.
This paper proposes a regularized GCov estimator to address the invertibility issue of the variance matrix in the objective function of the GCov estimator due to the existence of many variables or consideration of many nonlinear transformations. We consider the Ridge-type regularization for the sample variance matrix, and we show that RGCov is consistent and has an asymptotically Normal distribution. Moreover, we provide conditions for the RGCov estimator to reach the semi-parametric efficiency bound.
We introduce the RGCov specification test on estimated residuals and the RNLSD test to test the absence of linear and nonlinear serial dependence in time series. We show that both tests, based on the regularization parameter converging to a nonzero constant, have an asymptotically a mixture of chi-square distributions and, based on the regularization parameter converging to zero, have an asymptotically a chi-square distribution with known degrees of freedom.
We explore the finite sample properties of the RGCov estimator based on vast simulation studies. We consider many variables and transformations that could cause an invertibility issue of the variance matrix. In both cases, the RGCov estimator performs significantly better than the GCov and diagonal GCov estimators. For an empirical illustration, we use the RGCov estimator to fit a mixed VAR model with roots inside and outside of the unit circle to 12 stock price series of green energy companies included in the Renixx index. Consequently, based on causal and noncausal components of mixed VAR(1), we proposed portfolios that perform better than the Renixx index in terms of cumulative rate of return.