Modelling Volatile Time Series with V-Transforms and Copulas

: An approach to the modelling of volatile time series using a class of uniformity-preserving transforms for uniform random variables is proposed. V-transforms describe the relationship between quantiles of the stationary distribution of the time series and quantiles of the distribution of a predictable volatility proxy variable. They can be represented as copulas and permit the formulation and estimation of models that combine arbitrary marginal distributions with copula processes for the dynamics of the volatility proxy. The idea is illustrated using a Gaussian ARMA copula process and the resulting model is shown to replicate many of the stylized facts of ﬁnancial return series and to facilitate the calculation of marginal and conditional characteristics of the model including quantile measures of risk. Estimation is carried out by adapting the exact maximum likelihood approach to the estimation of ARMA processes, and the model is shown to be competitive with standard GARCH in an empirical application to Bitcoin return data. for models using v-transforms and copula processes and show that they rival and often outperform forecast models from the extended GARCH family.


Introduction
In this paper, we show that a class of uniformity-preserving transformations for uniform random variables can facilitate the application of copula modelling to time series exhibiting the serial dependence characteristics that are typical of volatile financial return data. Our main aims are twofold: to establish the fundamental properties of v-transforms and show that they are a natural fit to the volatility modelling problem; to develop a class of processes using the implied copula process of a Gaussian ARMA model that can serve as an archetype for copula models using v-transforms. Although the existing literature on volatility modelling in econometrics is vast, the models we propose have some attractive features. In particular, as copula-based models, they allow the separation of marginal and serial dependence behaviour in the construction and estimation of models.
A distinction is commonly made between genuine stochastic volatility models, as investigated by Taylor (1994) and Andersen (1994), and GARCH-type models as developed in a long series of papers by Engle (1982), Bollerslev (1986), Ding et al. (1993), Glosten et al. (1993) and Bollerslev et al. (1994), among others. In the former an unobservable process describes the volatility at any time point while in the latter volatility is modelled as a function of observable information describing the past behaviour of the process; see also the review articles by Shephard (1996) and Andersen and Benzoni (2009). The generalized autoregressive score (GAS) models of Creal et al. (2013) generalize the observation-driven approach of GARCH models by using the score function of the conditional density to model time variation in key parameters of the time series model. The models of this paper have more in common with the observation-driven approach of GARCH and GAS but have some important differences.
In GARCH-type models, the marginal distribution of a stationary process is inextricably linked to the dynamics of the process as well as the conditional or innovation distribution; in most cases, it has no simple closed form. For example, the standard GARCH mechanism serves to create power-law behaviour in the marginal distribution, even when the innovations come from a lighter-tailed distribution such as Gaussian (Mikosch and Stȃricȃ 2000). While such models work well for many return series, they may not be sufficiently flexible to describe all possible combinations of marginal and serial dependence behaviour encountered in applications. In the empirical example of this paper, which relates to log-returns on the Bitcoin price series, the data appear to favour a marginal distribution with sub-exponential tails that are lighter than power tails and this cannot be well captured by standard GARCH models. Moreover, in contrast to much of the GARCH literature, the models we propose make no assumptions about the existence of second-order moments and could also be applied to very heavy-tailed situations where variance-based methods fail.
Let X 1 , . . . , X n be a time series of financial returns sampled at (say) daily frequency and assume that these are modelled by a strictly stationary stochastic process (X t ) with marginal distribution function (cdf) F X . To match the stylized facts of financial return data described, for example, by Campbell et al. (1997) and Cont (2001), it is generally agreed that (X t ) should have limited serial correlation, but the squared or absolute processes (X 2 t ) and (|X t |) should have significant and persistent positive serial correlation to describe the effects of volatility clustering.
In this paper, we refer to transformed series like (|X t |), in which volatility is revealed through serial correlation, as volatility proxy series. More generally, a volatility proxy series (T(X t )) is obtained by applying a transformation T : R → R which (i) depends on a change point µ T that may be zero, (ii) is increasing in X t − µ T for X t ≥ µ T and (iii) is increasing in µ T − X t for X t ≤ µ T .
Our approach in this paper is to model the probability-integral transform (PIT) series (V t ) of a volatility proxy series. This is defined by V t = F T(X) (T(X t )) for all t, where F T(X) denotes the cdf of T(X t ). If (U t ) is the PIT series of the original process (X t ), defined by U t = F X (X t ) for all t, then a v-transform is a function describing the relationship between the terms of (V t ) and the terms of (U t ). Equivalently, a v-transform describes the relationship between quantiles of the distribution of X t and the distribution of the volatility proxy T(X t ). Alternatively, it characterizes the dependence structure or copula of the pair of variables (X t , T(X t )). In this paper, we show how to derive flexible, parametric families of v-transforms for practical modelling purposes.
To gain insight into the typical form of a v-transform, let x 1 , . . . , x n represent the realized data values and let u 1 , . . . , u n and v 1 , . . . , v n be the samples obtained by applying the transformations v t = F (|X|) n (|x t |) and u t = F (X) n (x t ), where F (X) n (x) = 1 n+1 ∑ n t=1 I {x t ≤x} and F (|X|) n (x) = 1 n+1 ∑ n t=1 I {|x t |≤x} denote scaled versions of the empirical distribution functions of the x t and |x t | samples, respectively. The graph of (u t , v t ) gives an empirical estimate of the v-transform for the random variables (X t , |X t |). In the left-hand plot of Figure 1 we show the relationship for a sample of n = 1043 daily log-returns of the Bitcoin price series for the years 2016-2019. Note how the empirical v-transform takes the form of a slightly asymmetric 'V'.
The right-hand plot of Figure 1 shows the sample autocorrelation function (acf) of the data given by z t = Φ −1 (v t ) where Φ is the standard normal cdf. This reveals a persistent pattern of positive serial correlation which can be modelled by the implied ARMA copula. This pattern is not evident in the acf of the raw x t data in the centre plot.
To construct a volatility model for (X t ) using v-transforms, we need to specify a process for (V t ). In principle, any model for a series of serially dependent uniform variables can be applied to (V t ). In this paper, we illustrate concepts using the Gaussian copula model implied by the standard ARMA dependence structure. This model is particularly tractable and allows us to derive model properties and fit models to data relatively easily.
There is a large literature on copula models for time series; see, for example, the review papers by Patton (2012) and Fan and Patton (2014). While the main focus of this literature has been on cross-sectional dependencies between series, there is a growing literature on models of serial dependence. First-order Markov copula models have been investigated by Chen and Fan (2006), Chen et al. (2009) and Domma et al. (2009) while higher-order Markov copula models using D-vines are applied by Smith et al. (2010). These models are based on the pair-copula apporoach developed in Joe (1996), Cooke (2001, 2002) and Aas et al. (2009). However, the standard bivariate copulas that enter these models are not generally effective at describing the typical serial dependencies created by stochastic volatility, as observed by Loaiza-Maya et al. (2018).
and F (|X|) n denote versions of the empirical distribution function of the x t and |x t | values, respectively. The sample size is n = 1043 and the data are daily log-returns of the Bitcoin price for the years 2016-2019.
The paper is structured as follows. In Section 2, we provide motivation for the paper by constructing a symmetric model using the simplest example of a v-transform. The general theory of v-transforms is developed in Section 3 and is used to construct the class of VT-ARMA processes and analyse their properties in Section 4. Section 5 treats estimation and statistical inference for VT-ARMA processes and provides an example of their application to the Bitcoin return data; Section 6 presents the conclusions. Proofs may be found in the Appendix A.

A Motivating Model
Given a probability space (Ω, F , P), we construct a symmetric, strictly stationary process (X t ) t∈N\{0} such that, under the even transformation T(x) = |x|, the serial dependence in the volatility proxy series (T(X t )) is of ARMA type. We assume that the marginal cdf F X of (X t ) is absolutely continuous and the density f X satisfies f X (x) = f X (−x) for all x > 0. Since F X and F |X| are both continuous, the properties of the probability-integral (PIT) transform imply that the series (U t ) and (V t ) given by U t = F X (X t ) and V t = F |X| (|X t |) both have standard uniform marginal distributions. Henceforth, we refer to (V t ) as the volatility PIT process and (U t ) as the series PIT process.
Any other volatility proxy series that can be obtained by a continuous and strictly increasing transformation of the terms of (|X t |), such as (X 2 t ), yields exactly the same volatility PIT process. For example, ifṼ t = F X 2 (X 2 t ), then it follows from the fact that In this sense, we can think of classes of equivalent volatility proxies, such as (|X t |), (X 2 t ), (exp |X t |) and (ln(1 + |X t |)). In fact, (V t ) is itself an equivalent volatility proxy to (|X t |) since F |X| is a continuous and strictly increasing transformation.
The symmetry of f X implies that which implies that the relationship between the volatility PIT process (V t ) and the series PIT process (U t ) is given by where V (u) = |2u − 1| is a perfectly symmetric v-shaped function that maps values of U t close to 0 or 1 to values of V t close to 1, and values close to 0.5 to values close to 0. V is the canonical example of a v-transform. It is related to the so-called tent-map transformation T (u) = 2 min(u, 1 − u) by V (u) = 1 − T (u). Given (V t ), let the process (Z t ) be defined by setting Z t = Φ −1 (V t ) so that we have the following chain of transformations: We refer to (Z t ) as a normalized volatility proxy series. Our aim is to construct a process (X t ) such that, under the chain of transformations in (2), we obtain a Gaussian ARMA process (Z t ) with mean zero and variance one. We do this by working back through the chain.
The transformation V is not an injection and, for any V t > 0, there are two possible inverse values, 1 2 (1 − V t ) and 1 2 (1 + V t ). However, by randomly choosing between these values, we can 'stochastically invert' V to construct a random variable U t such that V (U t ) = V t , This is summarized in Lemma 1, which is a special case of a more general result in Proposition 4. Lemma 1. Let V be a standard uniform variable. If V = 0, set U = 1 2 . Otherwise, let U = 1 2 (1 − V) with probability 0.5 and U = 1 2 (1 + V) with probability 0.5. Then, U is uniformly distributed and V (U) = V.
This simple result suggests Algorithm 1 for constructing a process (X t ) with symmetric marginal density f X such that the corresponding normalized volatility proxy process (Z t ) under the absolute value transformation (or continuous and strictly increasing functions thereof) is an ARMA process. We describe the resulting model as a VT-ARMA process.
Algorithm 1: 1. Generate (Z t ) as a causal and invertible Gaussian ARMA process of order (p, q) with mean zero and variance one.

Form the volatility PIT process
3. Generate a process of iid Bernoulli variables (Y t ) such that P(Y t = 1) = 0.5.
4. Form the PIT process (U t ) using the transformation It is important to state that the use of the Gaussian process (Z t ) as the fundamental building block of the VT-ARMA process in Algorithm 1 has no effect on the marginal distribution of (X t ), which is F X as specified in the final step of the algorithm. The process (Z t ) is exploited only for its serial dependence structure, which is described by a family of finite-dimensional Gaussian copulas; this dependence structure is applied to the volatility proxy process. Figure 2 shows a symmetric VT-ARMA(1,1) process with ARMA parameters α 1 = 0.95 and β 1 = −0.85; such a model often works well for financial return data. Some intuition for this observation can be gained from the fact that the popular GARCH(1,1) model is known to have the structure of an ARMA(1,1) model for the squared data process; see, for example, McNeil et al. (2015) (Section 4.2) for more details.  Figure 2. Realizations of length n = 500 of (X t ) and (Z t ) for a VT-ARMA(1,1) process with a marginal Student t distribution with ν = 3 degrees of freedom and ARMA parameters α = 0.95 and β = −0.85. ACF plots for (X t ) and (|X t |) are also shown.

V-Transforms
To generalize the class of v-transforms, we admit two forms of asymmetry in the construction described in Section 2: we allow the density f X to be skewed; we introduce an asymmetric volatility proxy.
Definition 1 (Volatility proxy transformation and profile). Let T 1 and T 2 be strictly increasing, continuous and differentiable functions on is a volatility proxy transformation. The parameter µ T is the change point of T and the associated function g T : R + → R + , g T (x) = T −1 2 • T 1 (x) is the profile function of T.
By introducing µ T , we allow for the possibility that the natural change point may not be identical to zero. By introducing different functions T 1 and T 2 for returns on either side of the change point, we allow the possibility that one or other may contribute more to the volatility proxy. This has a similar economic motivation to the leverage effects in GARCH models (Ding et al. 1993); falls in equity prices increase a firm's leverage and increase the volatility of the share price.
Clearly, the profile function of a volatility proxy transformation is a strictly increasing, continuous and differentiable function on R + such that g T (x) = 0. In conjunction with µ T , the profile contains all the information about T that is relevant for constructing vtransforms. In the case of a volatility proxy transformation that is symmetric about µ T , the profile satisfies g T (x) = x.
The following result shows how v-transforms V = V (U) can be obtained by considering different continuous distributions F X and different volatility proxy transformations T of type (3). Proposition 1. Let X be a random variable with absolutely continuous and strictly increasing cdf F X on R and let T be a volatility proxy transformation. Let U = F X (X) and V = F T(X) (T(X)).
The result implies that any two volatility proxy transformations T andT which have the same change point µ T and profile function g T belong to an equivalence class with respect to the resulting v-transform. This generalizes the idea that T(x) = |x| and T(x) = x 2 give the same v-transform in the symmetric case of Section 2. Note also that the volatility proxy transformations T (V) and T (Z) defined by are in the same equivalence class as T since they share the same change point and profile function.
Definition 2 (v-transform and fulcrum). Any transformation V that can be obtained from Equation (4) by choosing an absolutely continuous and strictly increasing cdf F X on R and a volatility proxy transformation T is a v-transform. The value δ = F X (µ T ) is the fulcrum of the v-transform.

A Flexible Parametric Family
In this section, we derive a family of v-transforms using construction (4) by taking a tractable asymmetric model for F X using the construction proposed by Fernández and Steel (1998) and by setting µ T = 0 and g T (x) = kx ξ for k > 0 and ξ > 0. This profile function contains the identity profile g T (x) = x (corresponding to the symmetric volatility proxy transformation) as a special case, but allows cases where negative or positive returns contribute more to the volatility proxy. The choices we make may at first sight seem rather arbitrary, but the resulting family can in fact assume many of the shapes that are permissible for v-transforms, as we will argue.
Let f 0 be a density that is symmetric about the origin and let γ > 0 be a scalar parameter. Fernandez and Steel suggested the model This model is often used to obtain skewed normal and skewed Student distributions for use as innovation distributions in econometric models. A model with γ > 1 is skewed to the right while a model with γ < 1 is skewed to the left, as might be expected for asset returns. We consider the particular case of a Laplace or double exponential distribution f 0 (x) = 0.5 exp(−|x|) which leads to particularly tractable expressions.
It is remarkable that (7) is a uniformity-preserving transformation. If we set ξ = 1 and κ = 1, we get (8) is a very convenient special case, and we refer to it as the linear v-transform.
In Figure 3, we show the v-transform V δ,κ,ξ when δ = 0.55, κ = 1.4 and ξ = 0.65. We will use this particular v-transform to illustrate further properties of v-transforms and find a characterization.

Characterizing v-Transforms
It is easily verified that any v-transform obtained from (4) consists of two arms or branches, described by continuous and strictly monotonic functions; the left arm is decreasing and the right arm increasing. See Figure 3 for an illustration. At the fulcrum δ, we have V (δ) = 0. Every point u ∈ [0, 1] \ {δ} has a dual point u * on the opposite side of the fulcrum such that V (u * ) = V (u). Dual points can be interpreted as the quantile probability levels of the distribution of X that give rise to the same level of volatility.  (7). For any v-transform, if v = V (u) and u * is the dual of u, then the points (u, 0), (u, v), (u * , 0) and (u * , v) form the vertices of a square. For the given fulcrum δ, a v-transform can never enter the gray shaded area of the plot.
We collect these properties together in the following lemma and add one further important property that we refer to as the square property of a v-transform; this property places constraints on the shape that v-transforms can take and is illustrated in Figure 3.
V is continuous;
It is instructive to see why the square property must hold. Consider Figure 3 and fix a point u ∈ The properties in Lemma 2 could be taken as the basis of an alternative definition of a v-transform. In view of (9), it is clear that any mapping V that has these properties is a uniformity-preserving transformation. We can characterize the mappings V that have these properties as follows.
where Ψ is a continuous and strictly increasing distribution function on [0, 1].
Our arguments so far show that every v-transform must have the form (10). It remains to verify that every uniformity-preserving transformation of the form (10) can be obtained from construction (4), and this is the purpose of the final result of this section. This allows us to view Definition 2, Lemma 2, and the characterization (10) as three equivalent approaches to the definition of v-transforms.
Proposition 3. Let V be a uniformity-preserving transformation of the form (10) and F X a continuous distribution function. Then, V can be obtained from construction (4) using any volatility proxy transformation with change point µ T = F −1 X (δ) and profile Henceforth, we can view (10) as the general equation of a v-transform. Distribution functions Ψ on [0, 1] can be thought of as generators of v-transforms. Comparing (10) with (7), we see that our parametric family V δ,κ,ξ is generated by Ψ(x) = exp(−κ(−(ln x) ξ )). This is a 2-parameter distribution whose density can assume many different shapes on the unit interval including increasing, decreasing, unimodal, and bathtub-shaped forms. In this respect, it is quite similar to the beta distribution which would yield an alternative family of v-transforms. The uniform distribution function Ψ(x) = x gives the family of linear v-transforms V δ .
In applications, we construct models starting from the building blocks of a tractable v-transform V such as (7) and a distribution F X ; from these, we can always infer an implied profile function g T using (11). The alternative approach of starting from g T and F X and constructing V via (4) is also possible but can lead to v-transforms that are cumbersome and computationally expensive to evaluate if F X and its inverse do not have simple closed forms.

V-Transforms and Copulas
If two uniform random variables are linked by the v-transform V = V (U), then the joint distribution function of (U, V) is a special kind of copula. In this section, we derive the form of the copula, which facilitates the construction of stochastic processes using v-transforms.
To state the main result, we use the notation V −1 and V for the the inverse function and the gradient function of a v-transform V. Although there is no unique inverse V −1 (v) (except when v = 0), the fact that the two branches of a v-transform mutually determine each other allows us to define V −1 (v) to be the inverse of the left branch of the v-transform given by V −1 : Theorem 2. Let V and U be random variables related by the v-transform V = V (U).

1.
The joint distribution function of (U, V) is given by the copula 2.
Conditional on V = v, the distribution of U is given by 3.
. We note that this copula is related to a special case of the tent map copula family C T θ in Rémillard (2013) For the linear v-transform family, the conditional probability ∆(v) in (14) satisfies ∆(v) = δ. This implies that the value of V contains no information about whether U is likely to be below or above the fulcrum; the probability is always the same regardless of V. In general, this is not the case and the value of V does contain information about whether U is large or small.
Part (2) of Theorem 2 is the key to stochastically inverting a v-transform in the general case. Based on this result, we define the concept of stochastic inversion of a v-transform. We refer to the function ∆ as the conditional down probability of V.
is the stochastic inversion function of V.
The following proposition, which generalizes Lemma 1, allows us to construct general asymmetric processes that generalize the process of Algorithm 1.
Proposition 4. Let V and W be iid U(0, 1) variables and let V be a v-transform with stochastic inversion function V. If U = V −1 (V, W), then V (U) = V and U ∼ U(0, 1).
In Section 4, we apply v-transforms and their stochastic inverses to the terms of time series models. To understand the effect this has on the serial dependencies between random variables, we need to consider multivariate componentwise v-transforms of random vectors with uniform marginal distributions and these can also be represented in terms of copulas. We now give a result which forms the basis for the analysis of serial dependence properties. The first part of the result shows the relationship between copula densities under componentwise v-transforms. The second part shows the relationship under the componentwise stochastic inversion of a v-transform; in this case, we assume that the stochastic inversion of each term takes place independently given V so that all serial dependence comes from V .
Theorem 3. Let V be a v-transform and let U = (U 1 , . . . , U d ) and V = (V 1 , . . . , V d ) be vectors of uniform random variables with copula densities c U and c V , respectively.

1.
If where W 1 , . . . , W d are iid uniform random variables that are also independent of V 1 , . . . , V d , then

VT-ARMA Copula Models
In this section, we study some properties of the class of time series models obtained by the following algorithm, which generalizes Algorithm 1. The models obtained are described as VT-ARMA processes since they are stationary time series constructed using the fundamental building blocks of a v-transform V and an ARMA process.
We can add any marginal behaviour in the final step, and this allows for an infinitely rich choice. We can, for instance, even impose an infinite-variance or an infinite-mean distribution, such as the Cauchy distribution, and still obtain a strictly stationary process for (X t ). We make the following definitions.
Definition 4 (VT-ARMA and VT-ARMA copula process). Any stochastic process (X t ) that can be generated using Algorithm 2 by choosing an underlying ARMA process with mean zero and variance one, a v-transform V, and a continuous distribution function F X is a VT-ARMA process. The process (U t ) obtained at the penultimate step of the algorithm is a VT-ARMA copula process.

Algorithm 2:
1. Generate (Z t ) as a causal and invertible Gaussian ARMA process of order (p, q) with mean zero and variance one.
4. Form the series PIT process (U t ) by taking the stochastic inverses U t = V −1 (V t , W t ).
5. Form the process (X t ) by setting X t = F −1 X (U t ) for some continuous cdf F X . Figure 4 gives an example of a simulated process using Algorithm 2 and the vtransform V δ,κ,ξ in (7) with κ = 0.9 and MA parameter ξ = 1.1. The marginal distribution is a heavy-tailed skewed Student distribution of type (6) with degrees-of-freedom ν = 3 and skewness γ = 0.8, which gives rise to more large negative returns than large positive returns. The underlying time series model is an ARMA(1,1) model with AR parameter α = 0.95 and MA parameter β = −0.85. See the caption of the figure for full details of parameters.
In the remainder of this section, we concentrate on the properties of VT-ARMA copula processes (U t ) from which related properties of VT-ARMA processes (X t ) may be easily inferred.

Stationary Distribution
The VT-ARMA copula process (U t ) of Definition 4 is a strictly stationary process since the joint distribution of (U t 1 , . . . , U t k ) for any set of indices t 1 < · · · < t k is invariant under time shifts. This property follows easily from the strict stationarity of the underlying ARMA process (Z t ) according to the following result, which uses Theorem 3.
An expression for the joint density facilitates the calculation of a number of dependence measures for the bivariate marginal distribution of (U t , U t+k ). In the bivariate case, the correlation matrix of the underlying Gaussian copula C Ga P(t,t+k) contains a single offdiagonal value ρ(k) and we simply write C Ga ρ(k) . The Pearson correlation of (U t , U t+k ) is given by This value is also the value of the Spearman rank correlation ρ S (X t , X t+k ) for a VT-ARMA process (X t ) with copula process (U t ) (since the Spearman's rank correlation of a pair of continuous random variables is the Pearson correlation of their copula). The calculation of (18) typically requires numerical integration. However, in the special case of the linear v-transform V δ in (8), we can get a simpler expression as shown in the following result. Proposition 6. Let (U t ) be a VT-ARMA copula process satisfying the assumptions of Proposition 5 with linear v-transform V δ . Let (Z t ) denote the underlying Gaussian ARMA process. Then, For the symmetric v-transform V 0.5 , Equation (19) obviously yields a correlation of zero so that, in this case, the VT-ARMA copula process (U t ) is a white noise with an autocorrelation function that is zero, except at lag zero. However, even a very asymmetric model with δ = 0.4 or δ = 0.6 gives ρ(U t , U t+k ) = 0.04ρ S (Z t , Z t+k ) so that serial correlations tend to be very weak.
When we add a marginal distribution, the resulting process (X t ) has a different autocorrelation function to (U t ), but the same rank autocorrelation function. The symmetric model of Section 2 is a white noise process. General asymmetric processes (X t ) are not perfect white noise processes but have only very weak serial correlation.

Conditional Distribution
To derive the conditional distribution of a VT-ARMA copula process, we use the vector notation U t = (U 1 , . . . , U t ) and Z t = (Z 1 , . . . , Z t ) to denote the history of processes up to time point t and u t and z t for realizations. These vectors are related by the componentwise transformation Z t = Φ −1 (V (U t )). We assume that all processes have a time index set given by t ∈ {1, 2, . . .}.
Proposition 7. For t > 1, the conditional density f U t |U t−1 (u | u t−1 ) is given by where µ t = E(Z t | Z t−1 = Φ −1 (V (u t−1 ))) and σ is the standard deviation of the innovation process for the ARMA model followed by (Z t ).
When (Z t ) is iid white noise µ t = 0, σ = 1 and (20) reduce to the uniform density f U t |U t−1 (u | u t−1 ) = 1 as expected. In the case of the first-order Markov AR(1) model Z t = α 1 Z t−1 + t , the conditional mean of Z t is µ t = α 1 Φ −1 (V (u t−1 )) and σ 2 = 1 − α 2 1 . The conditional density (20) can be easily shown to simplify to ) denotes the copula density derived in Proposition 5. In this special case, the VT-ARMA model falls within the class of first-order Markov copula models considered by Chen and Fan (2006), although the copula is new.
If we add a marginal distribution F X to the VT-ARMA copula model to obtain a model for (X t ) and use similar notational conventions as above, the resulting VT-ARMA model has conditional density with f U t |U t−1 as in (20). An interesting property of the VT-ARMA process is that the conditional density (21) can have a pronounced bimodality for values of µ t in excess of zero that is in high volatility situations where the conditional mean of Z t is higher than the marginal mean value of zero; in low volatility situations, the conditional density appears more concentrated around zero. This phenomenon is illustrated in Figure 4. The bimodality in high volatility situations makes sense: in such cases, it is likely that the next return will be large in absolute value and relatively less likely that it will be close to zero. The conditional distribution function of (X t ) is F X t |X t−1 (x | x t−1 ) = F U t |U t−1 (F X (x) | F X (x t−1 )) and hence the ψ-quantile x ψ,t of F X t |X t−1 can be obtained by solving For ψ < 0.5, the negative of this value is often referred to as the conditional (1 − ψ)-VaR (value-at-risk) at time t in financial applications.

Statistical Inference
In the copula approach to dependence modelling, the copula is the object of central interest and marginal distributions are often of secondary importance. A number of different approaches to estimation are found in the literature. As before, let x 1 , . . . , x n represent realizations of variables X 1 , . . . , X n from the time series process (X t ).
The semi-parametric approach developed by Genest et al. (1995) is very widely used in copula inference and has been applied by Chen and Fan (2006) to first-order Markov copula models in the time series context. In this approach, the marginal distribution F X is first estimated non-parametrically using the scaled empirical distribution function F (X) n (see definition in Section 1) and the data are transformed onto the (0, 1) scale. This has the effect of creating pseudo-copula data u t = rank(x t )/(n + 1) where rank(x t ) denotes the rank of x t within the sample. The copula is fitted to the pseudo-copula data by maximum likelihood (ML).
As an alternative, the inference-functions-for-margins (IFM) approach of Joe (2015) could be applied. This is also a two-step method although in this case a parametric model F X is estimated under an iid assumption in the first step and the copula is fitted to the data u t = F X (x t ) in the second step.
The approach we adopt for our empirical example is to first use the semi-parametric approach to determine a reasonable copula process, then to estimate marginal parameters under an iid assumption, and finally to estimate all parameters jointly using the parameter estimates from the previous steps as starting values.
We concentrate on the mechanics of deriving maximum likelihood estimates (MLEs). The problem of establishing the asymptotic properties of the MLEs in our setting is a difficult one. It is similar to, but appears to be more technically challenging than, the problem of showing consistency and efficiency of MLEs for a Box-Cox-transformed Gaussian ARMA process, as discussed in Terasaka and Hosoya (2007). We are also working with a componentwise transformed ARMA process, although, in our case, the transformation (X t ) → (Z t ) is via the nonlinear, non-increasing volatility proxy transformation T (Z) (x) in (5), which is not differentiable at the change point µ T . We have, however, run extensive simulations which suggests good behaviour of the MLEs in large samples.

Maximum Likelihood Estimation of the VT-ARMA Copula Process
We first consider the estimation of the VT-ARMA copula process for a sample of data u 1 , . . . , u n . Let θ (V) and θ (A) denote the parameters of the v-transform and ARMA model, respectively. It follows from Theorem 3 (part 2) and Proposition 5 that the loglikelihood for the sample u 1 , . . . , u n is simply the log density of the Gaussian copula under componentwise inverse v-transformation. This is given by where the first term L * is the log-likelihood for an ARMA model with a standard N(0,1) marginal distribution. Both terms in the log-likelihood (23) are relatively straightforward to evaluate. The evaluation of the ARMA likelihood L * (θ (A) | z 1 , . . . , z n ) for parameters θ (A) and data z 1 , . . . , z n can be accomplished using the Kalman filter. However, it is important to note that the assumption that the data z 1 , . . . , z n are standard normal requires a bespoke implementation of the Kalman filter, since standard software always treats the error variance σ 2 as a free parameter in the ARMA model. In our case, we need to constrain σ 2 to be a function of the ARMA parameters so that var(Z t ) = 1. For example, in the case of an ARMA(1,1) model with AR parameter α 1 and MA parameter β 1 , this means that σ 2 = σ 2 (α 1 , β 1 ) = (1 − α 2 1 )/(1 + 2α 1 β 1 + β 2 1 ). The constraint on σ 2 must be incorporated into the state-space representation of the ARMA model.
Model validation tests for the VT-ARMA copula can be based on residuals where z t denotes the implied realization of the normalized volatility proxy variable and where an estimate µ t of the conditional mean µ t = E(Z t | Z t−1 = z t ) may be obtained as an output of the Kalman filter. The residuals should behave like an iid sample from a normal distribution. Using the estimated model, it is also possible to implement a likelihood-ratio (LR) test for the presence of stochastic volatility in the data. Under the null hypothesis that θ (A) = 0, the log-likelihood (23) is identically equal to zero. Thus, the size of the maximized log-likelihood L( θ (V) , θ (A) ; u 1 , . . . , u n ) provides a measure of the evidence for the presence of stochastic volatility.

Adding a Marginal Model
If F X and f X denote the cdf and density of the marginal model and the parameters are denoted θ (M) , then the full log-likelihood for the data x 1 , . . . , x n is simply where the first term is the log-likelihood for a sample of iid data from the marginal distribution F X and the second term is (23). When a marginal model is added, we can recover the implied form of the volatility proxy transformation using Proposition 3. If δ is the estimated fulcrum parameter of the v-transform, then the estimated change point is µ T = F −1 X ( δ; θ (M) ) and the implied profile function is Note that is is possible to force the change point to be zero in a joint estimation of marginal model and copula by imposing the constraint F X (0; θ (M) ) = δ on the fulcrum and marginal parameters during the optimization. However, in our experience, superior fits are obtained when these parameters are unconstrained.

Example
We analyse n = 1043 daily log-returns for the Bitcoin price series for the period 2016-2019; values are multiplied by 100. We first apply the semi-parametric approach of Genest et al. (1995) using the log-likelihood (23) which yields the results in Table 1. Different models are referred to by VT(n)-ARMA(p, q), where (p, q) refers to the ARMA model and n indexes the v-transform: 1 is the linear v-transform V δ in (8); 3 is the three-parameter transform V δ,κ,ξ in (7); 2 is the two-parameter v-transform given by V δ,κ := V δ,κ,1 . In unreported analyses, we also tried the three-parameter family based on the beta distribution, but this had negligible effect on the results.
The column marked L gives the value of the maximized log-likelihood. All values are large and positive showing strong evidence of stochastic volatility in all cases. The model VT(1)-ARMA(1,0) is a first-order Markov model with linear v-transform. The fit of this model is noticeably poorer than the others suggesting that Markov models are insufficient to capture the persistence of stochastic volatility in the data. The column marked SW contains the p-value for a Shapiro-Wilks test of normality applied to the residuals from the VT-ARMA copula model; the result is non-significant in all cases. Table 1. Analysis of daily Bitcoin return data 2016-2019. Parameter estimates, standard errors (below estimates) and information about the fit: SW denotes Shapiro-Wilks p-value; L is the maximized value of the log-likelihood and AIC is the Akaike information criterion. According to the AIC values, the VT(2)-ARMA(1,1) is the best model. We experimented with higher order ARMA processes, but this did not lead to further significant improvements. Figure 5 provides a visual of the fit of this model. The pictures in the panels show the QQplot of the residuals against normal, acf plots of the residuals and squared residuals and the estimated conditional mean process ( µ t ), which can be taken as an indicator of high and low volatility periods. The residuals and absolute residuals show very little evidence of serial correlation and the QQplot is relatively linear, suggesting that the ARMA filter has been successful in explaining much of the serial dependence structure of the normalized volatility proxy process. We now add various marginal distributions to the VT(2)-ARMA(1,1) copula model and estimate all parameters of the model jointly. We have experimented with a number of location-scale families including Student-t, Laplace (double exponential), and a double-Weibull family which generalizes the Laplace distribution and is constructed by taking back-to-back Weibull distributions. Estimation results are presented for these three distributions in Table 2. All three marginal distributions are symmetric around their location parameters µ, and no improvement is obtained by adding skewness using the construction of Fernández and Steel (1998) described in Section 3.1; in fact, the Bitcoin returns in this time period show a remarkable degree of symmetry. In the table, the shape and scale parameters of the distributions are denoted η and σ, respectively; in the case of Student, an infinite-variance distribution with degree-of-freedom parameter η = 1.94 is fitted, but this model is inferior to the models with Laplace and double-Weibull margins; the latter is the favoured model on the basis of AIC values. Figure 6 shows some aspects of the joint fit for the fully parametric VT(2)-ARMA(1,1) model with double-Weibull margin. A QQplot of the data against the fitted marginal distribution confirms that the double-Weibull is a good marginal model for these data. Although this distribution is sub-exponential (heavier-tailed than exponential), its tails do not follow a power law and it is in the maximum domain of attraction of the Gumbel distribution (see, for example, McNeil et al. 2015, Chapter 5). Table 2. VT(2)-ARMA(1,1) model with three different margins: Student-t, Laplace, double Weibull. Parameter estimates, standard errors (alongside estimates) and information about the fit: SW denotes Shapiro-Wilks p-value; L is the maximized value of the log-likelihood and AIC is the Akaike information criterion.

Student
Laplace Using (26), the implied volatility proxy profile function g T can be constructed and is found to lie just below the line y = x as shown in the upper-right panel. The change point is estimated to be µ T = 0.06. We can also estimate an implied volatility proxy transformation in the equivalence class defined by g T and µ T . We estimate the transformation (M) ))). In the lower-left panel of Figure 6, we show the empirical v-transform formed from the data (x t , T(x t )) together with the fitted parametric v-transform V θ (V) . We recall from Section 1 that the empirical v-transform is the plot (u t , v t ) where u t = F (X) n (x t ) and v t = F ( T(X)) n ( T(x t )). The empirical v-transform and the fitted parametric v-transform show a good degree of correspondence. The lowerright panel of Figure 6 shows the volatility proxy transformation T(x) as a function of x superimposed on the points (x t , Φ −1 (v t )). Using the curve, we can compare the effects of, for example, a log-return (× 100) of -10 and a log-return of 10. For the fitted model, these are 1.55 and 1.66 showing that the up movement is associated with slightly higher volatility.
As a comparison to the VT-ARMA model, we fit standard GARCH(1,1) models using Student-t and generalized error distributions for the innovations; these are standard choices available in the popular rugarch package in R. The generalized error distribution (GED) contains normal and Laplace as special cases as well as a model that has a similar tail behaviour to Weibull; note, however, that, by the theory of Mikosch and Stȃricȃ (2000), the tails of the marginal distribution of the GARCH decay according to a power law in both cases. The results in Table 3 show that the VT(2)-ARMA(1,1) models with Laplace and double-Weibull marginal distributions outperform both GARCH models in terms of AIC values.  Figure 7 shows the in-sample 95% conditional value-at-risk (VaR) estimate based on the VT(2)-ARMA(1,1) model which has been calculated using (22). For comparison, a dashed line shows the corresponding estimate for the GARCH(1,1) model with GED innovations.
Finally, we carry out an out-of-sample comparison of conditional VaR estimates using the same two models. In this analysis, the models are estimated daily throughout the 2016-2019 period using a 1000-day moving data window and one-step-ahead VaR forecasts are calculated. The VT-ARMA model gives 47 exceptions of the 95% VaR and 11 exceptions of the 99% VaR, compared with expected numbers of 52 and 10 for a 1043 day sample, while the GARCH model leads to 57 and 12 exceptions; both models pass binomial tests for these exception counts. In a follow-up paper (Bladt and McNeil 2020), we conduct more extensive out-of-sample backtests for models using v-transforms and copula processes and show that they rival and often outperform forecast models from the extended GARCH family.

Conclusions
This paper has proposed a new approach to volatile financial time series in which v-transforms are used to describe the relationship between quantiles of the return distribution and quantiles of the distribution of a predictable volatility proxy variable. We have characterized v-transforms mathematically and shown that the stochastic inverse of a v-transform may be used to construct stationary models for return series where arbitrary marginal distributions may be coupled with dynamic copula models for the serial dependence in the volatility proxy.
The construction was illustrated using the serial dependence model implied by a Gaussian ARMA process. The resulting class of VT-ARMA processes is able to capture the important features of financial return series including near-zero serial correlation (white noise behaviour) and volatility clustering. Moreover, the models are relatively straightforward to estimate building on the classical maximum-likelihood estimation of an ARMA model using the Kalman filter. This can be accomplished in the stepwise manner that is typical in copula modelling or through joint modelling of the marginal and copula process. The resulting models yield insights into the way that volatility responds to returns of different magnitude and sign and can give estimates of unconditional and conditional quantiles (VaR) for practical risk measurement purposes.
There are many possible uses for VT-ARMA copula processes. Because we have complete control over the marginal distribution, they are very natural candidates for the innovation distribution in other time series models. For example, they could be applied to the innovations of an ARMA model to obtain ARMA models with VT-ARMA errors; this might be particularly appropriate for longer interval returns, such as weekly or monthly returns, where some serial dependence is likely to be present in the raw return data.
Clearly, we could use other copula processes for the volatility PIT process (V t ). The VT-ARMA copula process has some limitations: the radial symmetry of the underlying Gaussian copula means that the serial dependence between large values of the volatility proxy must mirror the serial dependence between small values; moreover, this copula does not admit tail dependence in either tail and it seems plausible that very large values of the volatility proxy might have a tendency to occur in succession.
To extend the class of models based on v-transforms, we can look for models for the volatility PIT process (V t ) with higher dimensional marginal distributions given by asymmetric copulas with upper tail dependence. First-order Markov copula models as developed in Chen and Fan (2006) can give asymmetry and tail dependence, but they cannot model the dependencies at longer lags that we find in empirical data. D-vine copula models can model higher-order Markov dependencies and Bladt and McNeil (2020) show that this is a promising alternative specification for the volatility PIT process.
Funding: This research received no external funding.

Data Availability Statement:
The analyses were carried out using R 4.0.2 (R Core Team, 2020) and the tscopula package (Alexander J. McNeil and Martin Bladt, 2020) available at https://github.com/ ajmcneil/tscopula. The full reproducible code and the data are available at https://github.com/ ajmcneil/vtarma.

Acknowledgments:
The author is grateful for valuable input from a number of researchers including Hansjoerg Albrecher, Martin Bladt, Valérie Chavez-Demoulin, Alexandra Dias, Christian Genest, Michael Gordy, Yen Hsiao Lok, Johanna Nešlehová, Andrew Patton, and Ruodu Wang. Particular thanks are due to Martin Bladt for providing the Bitcoin data and advice on the data analysis. The paper was completed while the author was a guest at the Forschungsinstitut für Mathematik (FIM) at ETH Zurich.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. Proofs
Appendix A.1. Proof of Proposition 1 We observe that, for x ≥ 0, {X t ≤ µ T } ⇐⇒ {U ≤ F X (µ T )} and in this case The cumulative distribution function F 0 (x) of the double exponential distribution is equal to 0.5e x for x ≤ 0 and 1 − 0.5e −x if x > 0. It is straightforward to verify that When g T (x) = kx ξ , we obtain for u ≤ δ that For u > δ, we make a similar calculation.
Appendix A.3. Proof of Theorem 1 It is easy to check that Equation (10) fulfills the list of properties in Lemma 2. We concentrate on showing that a function that has these properties must be of the form (10). It helps to consider the picture of a v-transform in Figure 3. Consider the lines v = 1 − u and v = δ − u for u ∈ [0, δ]. The areas above the former and below the latter are shaded gray.
The left branch of the v-transform must start at (0, 1), end at (δ, 0), and lie strictly between these lines in (0, δ). Suppose, on the contrary, that v = V (u) ≤ δ − u for u ∈ (0, δ). This would imply that the dual point u * given by u * = u + v satisfies u * ≤ δ which contradicts the requirement that u * must be on the opposite side of the fulcrum. Similarly, if v = V (u) ≥ 1 − u for u ∈ (0, δ), then u * ≥ 1 and this is also not possible; if u * = 1, then u = 0, which is a contradiction.
Thus, the curve that links (0, 1) and (δ, 0) must take the form where Ψ(0) = 0, Ψ(1) = 1 and 0 < Ψ(x) < 1 for x ∈ (0, 1). Clearly, Ψ must be continuous to satisfy the conditions of the v-transform. It must also be strictly increasing. If it were not, then the derivative would satisfy V (u) ≥ −1, which is not possible: if at any point u ∈ (0, δ), we have V (u) = −1, then the opposite branch of the v-transform would have to jump vertically at the dual point u * , contradicting continuity; if V (u) > −1, then V would have to be a decreasing function at u * , which is also a contradiction. Thus, Ψ fulfills the conditions of a continuous, strictly increasing distribution function on [0, 1], and we have established the necessary form for the left branch equation. To find the value of the right branch equation at u > δ, we invoke the square property. Since V (u) = V (u * ) = V (u − V (u)), we need to solve the equation x = V (u − x) for x ∈ [0, 1] using the formula for the left branch equation of V. Thus, we solve x = 1 − u + x − (1 − δ)Ψ( u−x δ ) for x, and this yields the right branch equation as asserted.

2.
We can write P(U ≤ u, V ≤ v) = C(u, v), where C is the copula given by (12). It follows from the basic properties of a copula that This is the distribution function of a binomial distribution, and it must be the case that ∆(v) = − d dv V −1 (v). Equation (14) follows by differentiating the inverse. 3.
The value of Spearman's rho ρ S (Z t , Z t+k ) for the bivariate Gaussian distribution is well known; see, for example, McNeil et al. (2015).
The Gaussian copula density is given in general by where Z is a multivariate Gaussian vector with standard normal margins and correlation matrix P. Hence, it follows that we can write where f Z t |Z t−1 is the conditional density of the ARMA process, from which (20) follows easily.