{'fraction_non_alphanumeric': 0.04822745592630165, 'fraction_numerical': 0.015902174211060237, 'mean_word_length': 4.3363569861009506, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulation of Multivariate Non-Gaussian Autoregressive Time Series with Given Autocovariance and Marginals

A semi-analytic method is proposed for the generation of realizations of a multivariate process of a given linear correlation structure and marginal distribution. This is an extension of a similar method for univariate processes, transforming the autocorrelation of the non-Gaussian process to that of a Gaussian process based on a piece-wise linear marginal transform from non-Gaussian to Gaussian marginal. The extension to multivariate processes involves the derivation of the autocorrelation matrix from the marginal transforms, which determines the generating vector autoregressive process. The effectiveness of the approach is demonstrated on systems designed under different scenarios of autocovariance and marginals.


Introduction
Many real time series cannot be considered to be Gaussian and do not fit in the framework of standard linear analysis. A related problem is the generation of non-Gaussian time series with given linear correlation structure and marginal distribution. This problem arises mainly in stochastic simulation and randomization testing.
In stochastic simulation and particularly in the framework of input modelling, the problem occurs when there are dependencies among random variables that constitute the inputs to the simulation model. For time series data, this leads to the simulation of univariate or multivariate stochastic linear processes of given marginal distributions and lagged correlation structure [1,2,3,4]. This setting is met in many applications, ranging from manufacturing systems [5], to medical treatment [6], internet traffic [7], bird flocking [8], floods [9], and ocean temperature [10].
Randomization testing has been used to investigate nonlinear dependencies in the time series, where the null hypothesis is that the underlying process is linear stochastic and the test statistic is nonlinear. For many nonlinear statistics the null distribution is not known and a common approach, known as surrogate data test, is to form the empirical null distribution from the values of the test statistic computed on randomized time series consistent to the null hypothesis. The latter requires that the surrogate time series preserve the original marginal distribution and linear correlation structure [11,12,13,14]. The surrogate data test for nonlinearity has been mainly developed for univariate time series and has been applied in many fields for the investigation of nonlinear dynamics and chaos, such as finance [15,16,17], geophysics [18,19,20] and physiology [21,22]. For multivariate time series, there are few approaches approximating the linear correlation structure in the frequency domain using the cross-power spectrum [23,24,25].
Though the problem is the same, solutions were proposed independently in the two areas of stochastic stimulation and surrogate data testing. In stochastic simulation, the problem was postulated as generating time series of arbitrary length with a given autocorrelation and marginal distribution. In a series of works, the algorithm for univariate time series called autoregressive-to-anything (ARTA) was developed, modified and tested [26,1,27,28], and was further extended for multivariate time series, known as Vector ARTA (VARTA) [29]. The method relies on solving numerically a double integral expression for the transform of the product moment of two Gaussian variables to that of variables of arbitrary marginals, the latter being the lagged variables of the linear stochastic process of arbitrary marginal distribution. The computations may be simplified approximating the marginal distribution with the Johnson translation system of distributions [29], and the generalized Pareto distribution [9].
For the surrogate data test for nonlinearity, independently of the ARTA approaches, different algorithms were developed to generate time series that match the marginal distribution and correlation structure of the examined time series. All these methods use data randomization, matching exactly the sample marginal distribution, and approximating the linear correlation either in the frequency domain using the Fourier transform [11,12], refined further using wavelet transform [30], or in the time domain, correcting the autocorrelation function [31], or finding an appropriate Gaussian autoregressive process, called statically transformed autoregressive process (STAP) [14]. STAP and ARTA are similar in that both methods attempt to form the transform from Gaussian autocorrelation to the given autocorrelation, but ARTA uses numerical solution of the double integral form, while STAP uses parametric approximation, originally polynomial [14], and then piece-wise linear [32].
Here, we extend the piece-wise approximation in STAP to estimate the auto-and cross-correlation of the multivariate time series, and we term the method as vector STAP (VSTAP). We demonstrate the performance of VSTAP on different simulated multivariate stochastic processes. The structure of the paper is as follows. In Section 2, we give the background and briefly discuss STAP using piece-wise approximation for univariate time series, and in Section 3, we present VSTAP for multivariate time series. In Section 4, we show the results of simulations on different multivariate stochastic processes, and we conclude in Section 5.

Univariate Time Series with given Marginal Distribution and Autocorrelation
We start with the univariate case, and suppose it is given a univariate time series {x t } n t=1 with marginal distribution F X (x) and sample autocorrelation function r X (τ ), where τ is the time lag. Equivalently, instead of r X (τ ) the power spectrum S X (f ) may be considered, where f is the frequency. The problem is to generate a time series {x * t } N t=1 , where N may be different from n, fulfilling the following two conditions: for a sufficiently large P . The second condition for the preservation of the linear correlation structure can be equivalently given in terms of power spectrum, S X * (f ) = S X (f ) for all frequencies f . In stochastic simulation, the problem may be postulated without reference to a specific time series {x t } n t=1 but only to the given F X (x) and r X (τ ).

Proposed solutions
The solutions we consider here match exactly the first condition in eq.(1) and approximate the second condition in eq.(2). All solutions make use of the marginal transform from Gaussian to the given distribution and the inverse transform assuming a variable Z following the standard Gaussian distribution with cumulative density function (cdf) Φ. There are two main approaches for the solution: the constrained realization approach, where the objective is to transform a random time series in order to match the given two conditions, and the typical realization approach, attempting to find a generating process that fulfills the two conditions. Though there has been some evidence in favor of constrained realization for hypothesis testing, as is the case with the surrogate data test for nonlinearity [33], it requires a time series {x t } n t=1 is given, from which the marginals and correlation structure are derived, and it gives another time series {x * t } n t=1 of the same length. On the other hand, the typical realization approach gives more insight onto the underlying process and has also the advantage that it can take as input only the marginals and autocorrelations, and can generate time series of any length, e.g. see [34]. The amplitude adjusted Fourier transform (AAFT) [11] and the iterated AAFT (IAAFT) [12] transform a time series to approximate the power spectrum and are constrained realization approaches. On the other hand, ARTA and STAP attempt to identify the process that generates realizations possessing the given autocorrelation and are therefore typical realization approaches. The two latter approaches are decomposed in the same four steps: 1. Starting with the marginal transform in eq.(3), the transform ψ τ from the Gaussian autocorrelation r Z (τ ) to the given autocorrelation r X (τ ) is determined for each lag τ , r X (τ ) = ψ τ (r Z (τ )), and the solution for each r Z (τ ) is obtained. 2. For a given order P , the coefficients of an AR(P ) process are computed from the autocorrelations r Z (τ ), τ = 1, . . . , P , using the Yule-Walker equations [35, Sec. 7.1]. 3. A Gaussian time series {z * t } N t=1 of a given length N is generated by the AR(P ) process. 4. The Gaussian time series is transformed to obtain the given marginal distribution,

The method of statically transformed autoregressive process
ARTA and STAP differ only in the first step. In ARTA, the double integral form for ψ τ is derived based on the marginal transform [1] where φ(z t , z t−τ , r Z (τ )) is the bivariate standard Gaussian probability density function (pdf), andx and s 2 are the sample mean and variance of X, respectively. To solve eq.(5) with respect to r Z (τ ), the double integral form is solved numerically (for an enhanced numerical solution, see [36]). On the other hand, STAP uses a parametric approximation of ψ τ . The original STAP in [14] uses polynomial approximation of the marginal transform in eq.(3), resulting in a polynomial form for ψ τ . We found that linear piecewise approximation gives a better solution of r X (τ ) = ψ τ (r Z (τ )) with respect to r Z (τ ) [32]. The linear piece-wise approximation of the marginal transform is comprised of first degree polynomials at each of m segments defined by m − 1 breakpoints a k , k = 1, . . . , m − 1 The partition of Z and the linear piece-wise function in (6) determine a partition Then the product moment for a lag τ is where P = P (a k−1 , a k , a l−1 , a l ; ρ Z (τ )) is the probability of (Z t , Z t−τ ) being in the region [a k−1 , a k ] × [a l−1 , a l ], µ 1,0 and µ 0,1 are the first order marginal moments and µ 1,1 the product moment of the doubly truncated Gaussian variables . Substituting the expressions for the moments of the joint doubly truncated Gaussian distribution, we get an analytic form for ψ τ (for details, see [32]). We found that best results are obtained when the breakpoints divide the standard Gaussian domain into equiprobable intervals. We also found that the constraint of continuity on the linear piece-wise function does not affect substantially the approximation of the correlation transform. Therefore we independently estimate the coefficients of the linear function at each interval rather than using linear splines. It is noted that the possible violation of continuity does not affect the monotonicity of the piece-wise linear function, which is always maintained. In a comparative study in [37], it was shown that r X * (τ ) from STAP with the polynomial approximation estimates r X (τ ) without bias, as opposed to AAFT, IAAFT and a model bootstrap approach, but with a larger variance than IAAFT, which decreases with the increase of n. The linear piecewise approximation decreases further the variance of r X * (τ ) from STAP and simulations in [32] showed that it gives more accurate Gaussian correlation estimation, and therefore we adopt it in the extension of STAP for multivariate time series presented below.

The Method of Vector Statically Transformed Autoregressive Process
For K multivariate time series, the problem involves the marginal distributions of all variables X 1 , . . . , X K , and the lagged cross-correlation for all pairs (X i , X j ) in addition to the autocorrelations for each X i . Given the marginal distributions F Xi (x), i = 1, . . . , K, and lagged correlations r Xi,Xj (τ ), the problem is to find a multivariate time series {x * i,t } N t=1 , i = 1, . . . , K, fulfilling the two following conditions: for a sufficiently large P . Alternatively, considering the linear structure in the frequency domain the second condition can be postulated in terms of cross-power spectrum, for all frequencies f . The first condition on the marginal distributions does not really complicate the solution for a proper {x * i,t } N t=1 as the marginal transform can be applied separately for each of the K variables. On the other hand, the preservation of the linear correlation structure is far more difficult to achieve than for the univariate case because in addition to the autocorrelation r Xi,Xi (τ ) (denoted r X (τ ) in the univariate case) the lagged cross-correlations r Xi,Xj (τ ), i = j, have to be matched as well (the same holds for the cross-power spectrum).

The typical realization approach
The typical realization approach is similar to that for univariate time series and is decomposed in the following four steps: 1. The correlation transform for each pair of variables (X i , X j ) and lag τ is formed as r Xi,Xj (τ ) = ψ i,j,τ (r Zi,Zj (τ )) from the marginal transforms, and the solution for each r Zi,Zj (τ ) is obtained. 2. Given the lagged cross-and auto-correlations r Zi,Zj (τ ), i, j = 1, . . . , K, τ = 0, . . . , P (for a given order P ), the coefficient matrices of a vector autoregressive process of order P on K variables, VAR K (P ), are computed using the multivariate generalization of the Yule-Walker equations [35,Sec. 16.5]. 3. The Gaussian multivariate time series {z * i,t } N t=1 , i = 1, . . . , K is generated by the VAR K (P ) process.

Each of the Gaussian time series {z
. Albeit the similarity of the four steps above to these for the univariate case in Sec. 2, there are additional problems in their implementation, which will be discussed later in Sec. 3.3. Both the vector ARTA (VARTA) [29], and vector STAP (VSTAP), presented below, implement the four steps for any marginal distributions and cross-correlation structure. Further, both methods can be implemented to multivariate time series assuming only that they are continuous valued and stationary. The main difference in VSTAP and VARTA is that in step 1, VARTA derives r Zi,Zj (τ ) solving numerically the double integral form as in eq.(5), while VSTAP approximates the marginal transform from Gaussian to the given sample distribution with a linear piece-wise function, given in eq.(6).

Implementation of VSTAP
VSTAP first fits the linear piece-wise function to the sample marginal transform, as for the univariate case (see eq. (6)). The breakpoints of the piece-wise function divide the standard Gaussian domain to equiprobable intervals and are thus the same for all variables, and only the coefficients of the linear piece-wise function are different for each variable Z i , denoted c i,k0 and c i,k1 , k = 1, . . . , m, for the constant term and the slope coefficient, respectively. The product moment for the pair of variables (X i , X j ) and lag τ is similar to that of eq.(7) for the univariate case , µ 1,0 and µ 0,1 are the first order marginal moments and µ 1,1 the product moment of the doubly truncated bivariate standard Gaussian distribution on [a k−1 , a k ] × [a l−1 , a l ]. The first order marginal moments in eq.(10) have the same expression for all i = 1, . . . , m, e.g. µ 1,0 is ϕ(u) du and ρ = ρ Zi,Zj (τ ). The product moment in eq.(10) is The probability P and the moments are functions of ρ Zi,Zj (τ ), so after substitution we get an analytic form forψ i,j,τ that approximates the true ψ i,j,τ based on linear piece-wise marginal transforms [32]. The functionψ i,j,τ is invertible, given that the piece-wise approximation in eq. (6) is also invertible [29, Theorem 3.4] (the theorem is based on the monotonicity of the cdf F Xi , and thus it applies to the piece-wise function fitted to F Xi being monotonic). However, the expression ofψ i,j,τ is very complicated and the closed form solution for its inverse cannot be obtained. We therefore use an iterative process to obtain the solution r Zi,Zj (τ ), approximating ρ Zi,Zj (τ ), so thatψ i,j,τ (r Zi,Zj (τ )) matches r Xi,Xj (τ ) at an arbitrary accuracy ǫ.
The iterative process needs a starting value for r Zi,Zj (τ ), and an appropriate value is given by the so-called naive correlation coefficient of X i,t and X j,t−τ . This is the Pearson correlation coefficient of the marginal transforms to Gaussian of X i,t and X j,t−τ according to eq.(4), r Zi, [38]. The main steps of the iterative algorithm are the following: 1. Begin with the naive correlation coefficient of X i,t and X j,t−τ as the starting value r Zi,Zj (τ ) (0) .

At each iteration
Otherwise the input for the next step is updated as r Zi,Zj (τ ) (h+1) = r Zi,Zj (τ ) (h) + δr Xi,Xj (τ ) and the computations are repeated from step 2.
The iterative algorithm is of the type of simple fixed point iteration, and thus has a linear rate of convergence. It requires only one starting value as opposed to bracketing methods (e.g. false position, secant, and bisection) requiring two starting values close and at each side of the solution [39,Chp.2]. The starting value of the naive correlation coefficient is close to the desired solution [38], and actually coincides with it in the case (X i,t , X j,t−τ ) is obtained from monotonic marginal transforms of a bivariate Gaussian variable pair [32]. At each iteration, we set the increment in the Gaussian correlation equal to the deviation in the target correlation δr Xi,Xj (τ ) because the Gaussian correlation and the target correlation are at the same amplitude level, and moreover it holds r Xi ,Xj (τ ) ≤ |r Zi,Zj (τ )| [40, p.600]. The linear convergence is guaranteed by the monotonicity of ψ i,j,τ and that r Xi,Xj (τ ) ≤ |r Zi,Zj (τ )|. In practice, we have found that the monotonicity ofψ i,j,τ may not hold at the edges of the interval [−1, 1]. For this, when |r Zi,Zj (τ ) (h) | is larger than a threshold close to one (we set the threshold to 0.9 to be on the safe side for all practical purposes) the condition of monotonicity is checked and if it is not satisfied, binary search is applied in the interval formed by the threshold and the current value. The number of iterations (including the binary search) depends on the given accuracy ǫ, but in any case the closed form expression ofψ i,j,τ (r Zi,Zj (τ )) makes the algorithm very time effective.

Implementation complications
A known problem with any correlation transform is that there may not be a feasible solution for a particular given correlation and marginals [26,41]. The domain of r Xi,Xj (τ ) for which a solution r Zi,Zj (τ ) ∈ [−1, 1] can be obtained is a subset of [−1, 1]. For a bivariate sample {x t , y t } n t=1 of (X, Y ) (in our case X = X i,t and Y = X j,t−τ ), the subset is formed by the minimum feasible correlation r and the maximum feasible correlation r where {x (t) } n t=1 is the ordered sample of X,x and s X are the sample mean and standard deviation [42]. For correlation matrices it is more difficult to determine the condition for feasibility, but there are computational procedures to check whether a given correlation matrix is feasible [41]. In the case of multivariate time series, the check for feasibility extends to the correlation matrix of any vector variable with components from the set {X i,t−τ | i = 1, . . . , K, τ = 0, . . . , P }. Even if all the correlation matrices are feasible, the solution for the corresponding Gaussian correlation matrices may not be valid, i.e. the matrices may not all be positive semi-definite. The problem arises because the components of each Gaussian correlation matrix are computed independently.
Instead of checking positive semidefiniteness for each possible Gaussian correlation matrix, the validity of the derived Gaussian lagged correlations can be tested collectively by checking for positive semidefiniteness of the full Gaussian correlation matrix where R Z (τ ) is the lagged correlation matrix for lag τ , τ = 0, . . . , P (P is the maximum lag) For correlation matrices, there are techniques to modify the matrix so as to be positive semidefinite, such as replacing negative eigenvalues with zero, or better with a slightly positive value in order to make it positive definite [43,41]. However, these techniques cannot be applied directly to the correlation matrix R Z in eq.(11) because it is comprised of repeated blocks being the lagged correlation matrices in eq.(12). Thus a change in an eigenvalue of R Z will alter its structure and make previously identical blocks differ. Our solution to this problem is to introduce a two-stage iterative procedure, where in the first stage we render positive definiteness of R Z and in the second stage we regain the structure in eq. (11). At each iteration, in the first stage we set the negative or zero eigenvalues of R Z to a slightly positive value and obtain a positive definite matrix but with altered components. In the second stage, for each repeated component in the form in eq.(11), we take the average of the values at the corresponding entries. For example, r Z1,Z2 (1) occurs in the entry (1, 2) of R Z (1) and entry (2, 1) of R Z (−1), so it occurs in all blocks R Z (1) and R Z (−1) of R Z , i.e. 2(P − 1) entries in total. Replacing all the repeated entries with the same average value gains back the correct structure of R Z but may cause R Z not to be positive definite, and the same two-stage procedure is then repeated. We have not worked out a proof for the convergence of this iterative procedure, but we found that a positive definite matrix R Z of the form in eq.(11) could be obtained after few steps.
According to step 2 of the typical realization approach in Sec. 3.1, given R Z the coefficient matrices A 1 , . . . , A P of the VAR K (P ) Gaussian process {Z 1,t , . . . , Z K,t } are computed from the multivariate generalization of the Yule-Walker equations [35,Sec. 16.5], where VAR K (P ) is expressed as and e t is uncorrelated process, here assumed to be Gaussian with unit covariance matrix. The positive definiteness of R Z in eq.(11) is important because it determines the stationarity of the VAR K (P ) Gaussian process {Z 1,t , . . . , Z K,t }.
In terms of A 1 , . . . , A P , VAR K (P ) is stationary when the roots of the reverse characteristic polynomial |I K − A 1 z − A 2 z 2 − · · · − A P z P | = 0 lie outside the unit circle in the complex plane, or equivalently the eigenvalues of are all smaller than one in modulus, where I K is the unit matrix. This condition is fulfilled by the condition of positive semidefiniteness of the full correlation matrix of R Z in eq.(11). Having a stationary VAR K (P ) Gaussian process, we generate a stationary Gaussian multivariate time series {z * i,t } N t=1 , i = 1, . . . , K from eq.(13). To transform the Gaussian {z * i,t } N t=1 to possess the given marginals in step 4 of the typical realization approach in Sec. 3.1, we can use either the marginal transform in eq.(3) or the linear piece-wise approximation. The former gives {x * i,t } N t=1 with exactly the same marginals and possibly some inaccuracy in the lagged correlations inherited by the inaccuracy of the linear piece-wise fit, while the latter loses some accuracy in matching the marginals but gains more accuracy in matching the lagged cross-correlations. The choice depends on the application. For example, for the randomization test for nonlinearity (where also we have N = n), we would choose the first approach to assure that the randomized time series contain exactly the same values as the original time series. We apply the latter approach in our simulations.

Simulations and Results
We show the efficiency of VSTAP in generating multivariate time series that match given non-Gaussian marginals and various correlation structures. In all simulations we set the number of breakpoints in the piecewise approximation to 20. A pilot study on smaller number of breakpoints showed that the decreased accuracy of the piece-wise linear fit does not affect much the solution for r Zi,Zj (τ ). So, for very small time series, a smaller number of breakpoints can also be used. On the other hand, our simulations on larger numbers of breakpoints showed insignificant improvement in the accuracy of the piece-wise linear fit and the estimated r Zi,Zj (τ ). For the accuracy tolerance, we use the absolute error rather than a relative error and set ǫ = 0.00001. For a generated realization of a process, VSTAP runs for the sample marginal distributions F Xi and the sample lagged correlations r Xi ,Xj (τ ). Using the approach that matches exactly the marginals in step 4 of Sec. 3.1, we establish that the marginals in {x * i,t } n t=1 coincide with those of the given time series {x i,t } n t=1 (we use n = N in the simulations). Thus only results for r Xi,Xj (τ ) are shown.
We consider different stationary VAR processes. We start with the Gaussian VAR 2 (2) generating process (K = 2, P = 2) where the input white noise vector e t = [e 1,t , e 2,t ] ′ is uncorrelated and has component variances σ 2 e1 = 0.09 and σ 2 e2 = 0.04 [44, p.17]. The observed time series {x t } n t=1 , where x t = [x 1,t , x 2,t ] ′ , has altered (non-Gaussian) marginals, given as x 1,t = s a 1,t , x 2,t = s a 2,t , and we set a = 3 and a = 2 to have a monotonic and a non-monotonic marginal transform, respectively. An example of realizations of the processes for a = 3 and a = 2 are shown in Fig. 1. For a = 3, the marginal distributions for both X 1 and X 2 have large kurtosis, while for a = 2 both marginal distributions are strongly right skewed. The match of r Xi,Xj (τ ) with VSTAP for the monotonic and non-monotonic marginal transform can be seen in Fig. 2. The 100 generated multivariate time (τ) Figure 2: The match of r X i ,X j (τ ), i, j = 1, 2, using VSTAP for monotonic (a = 3, in (a)) and non-monotonic (a = 2, in (b)) marginal transform of the Gaussian VAR 2 (2) generating process. The black line is for the given r X i ,X j (τ ) and the grey (cyan online) lines are for r X * i ,X * j (τ ) from 100 realizations of VSTAP (N = 1024). The two dashed black lines denote the 95% Fisher confidence intervals of r X i ,X j (τ ). In (a), the expected correlation from the cubic marginal transform is displayed by thick stippled line.
series {x * i,t } 1024 t=1 have sample lagged correlations r X * i ,X * j (τ ) that spread around the given r Xi,Xj (τ ), i.e. r Xi,Xj (τ ) is within the distribution of r X * i ,X * j (τ ). The spread of r X * i ,X * j (τ ) is at the level of the spread expected for sample Gaussian correlation of the same N , as indicated in Fig. 2 by the dashed black lines denoting the 95% Fisher confidence intervals of r Xi,Xj (τ ).
For the monotonic transform, it is possible to compare the given sample correlation and the VSTAP correlation with the theoretic correlation for the monotonically transformed Gaussian VAR 2 (2) process. For each Gaussian lagged correlation ρ Si,Sj (τ ) of the original Gaussian VAR 2 (2) obtained from the coefficients of VAR 2 (2) through the Yule-Walker equations, the transform x i,t = s 3 i,t , i = 1, 2, determines the correlation transform ρ Xi,Xj (τ ) = 2ρ Si,Sj (τ ) 3 /5 + 3ρ Si,Sj (τ )/5 [32]. It turns out that both the sample and VSTAP lagged correlations are close to the theoretic lagged correlations for (X 1 , X 1 ) and (X 1 , X 2 ) but differ for (X 2 , X 1 ) and (X 2 , X 2 ) (see Fig. 2a).
The mismatch of theoretical and sample correlation is more visible when the sample size N increases and the spread decreases, as shown in Fig. 3 for τ = 1. It is clearly shown in Fig. 3a that the mismatch when a = 3 occurs for ρ X2,X1 (1) (1) Figure 3: The match of r X i ,X j (1), i, j = 1, 2, against the time series length N (logarithmic scale) using VSTAP for monotonic (a = 3, in (a)) and non-monotonic (a = 2, in (b)) marginal transform of the Gaussian VAR 2 (2) generating process. The black line is for r X i ,X j (1) and the grey (cyan online) lines are for r X * i ,X * j (1) from 100 realizations of VSTAP. The dashed black lines denote the 95% Fisher confidence intervals of ρ X i ,X j (1). In (a), the expected correlation from the cubic marginal transform is displayed by thick stippled line. and ρ X2,X2 (1). For these cases, one can observe also difference in the sample correlations r X2,X1 (1) and r X2,X2 (1) and the VSTAP correlations r X * 2 ,X * 1 (1) and r X * 2 ,X * 2 (1), respectively, which is actually due to the bias in the estimation of the correlation matrix of the VAR process. Monte Carlo simulations on realizations of the Gaussian VAR 2 (2) process, without applying marginal transform and VSTAP, showed that the bias occurs even when the process is Gaussian. This bias is thus passed also to the VSTAP estimation of the given sample lagged correlations. For the monotonic marginal transform, we note that the mismatch of VSTAP for large N occurs only when there is bias, i.e. the sample correlation r Xi,Xj (1) differs from the theoretic correlation ρ Xi,Xj (1) (see Fig. 3a).
Finally, we report some results on the computational efficiency of VSTAP and compare it to VARTA. We focus on the computation of single components of the correlation matrix and consider a simplified setting of no auto-correlation, which is equivalent to having samples of correlated variables. In particular, we consider the example of a non-feasible Gaussian correlation matrix for three uniform variables with r X1,X2 (0) = −0.4, r X1,X3 (0) = 0.2 and r X2,X3 (0) = 0.8, first reported in [26]. The Gaussian correlation coefficients are given analytically from the expression r S1,S2 = 2 sin(πr X1,X2 /6) as r S1,S2 (0) = −0.4158, r S1,S3 (0)=0.2091 and r S2,S3 (0)=0.8135. The corresponding correlation matrix is not positive semidefinite and applying eigenvalue correction we derive the closest positive semidefinite correlation matrix with components r S1,S2 (0)=−0.4122, r S1,S3 (0) = 0.2062 and r S2,S3 (0) = 0.8065. Using these correlation coefficients, 1000 multivariate Gaussian samples are generated and transformed to uniform marginals applying the Gaussian cumulative density function. Since the marginal transform is monotonic we can evaluate VSTAP and VARTA (for this setting this is actually equivalent to the algorithm of normal to anything (NORTA) [36]) in matching r Xi,Xj (0) and r Si,Sj (0). The results for one of the three variable pairs is given in Figure 5, and similar are the results for the other two pairs. Both methods ran for the same accuracy level, given by a relative error of 0.001. Thus both methods match well r Xi,Xj (0), with r X * i ,X * j (0) of VARTA spreading evenly around r Xi,Xj (0), and r X * i ,X * j (0) of VSTAP being mostly concentrated at r Xi,Xj (0) and spread over at larger values, giving somewhat smaller root mean square error (RMSE). VSTAP gives also smaller RMSE in matching r Si,Sj (0), where both methods have much smaller accuracy (3) Figure 4: The match of r X i ,X j (τ ), i, j = 1, 2, against the time series length N (logarithmic scale) using VSTAP for monotonic (a = 3, first row) and non-monotonic (a = 2, second row) marginal transform of the Gaussian VAR 5 (4) generating process. The black line is for r X i ,X j (τ ) and the grey (cyan online) lines are for r X * i ,X * j (τ ) from 100 realizations of VSTAP, for the four largest r X i ,X j (τ ) found when N = 2 17 . The dashed black lines denote the 95% Fisher confidence intervals of r X i ,X j (τ ). In the first row the expected correlation from the cubic marginal transform is displayed by thick stippled line.
While both VSTAP and VARTA attain the same accuracy level in approximating r Xi,Xj (0), VSTAP succeeds this much faster 1 . As shown in Figure 6 for samples sizes n = 128, 256, 512, 1024, the computation time increases slowly with n for VSTAP and fast with VARTA. In particular, for VSTAP the scaling is logarithmic, and the fit for the computation time t in sec is t = −4.39 + 0.98 ln(n), while for VARTA t scales as a square of n, and the fit is t = −9.16n 1.98 . The  variance about the mean values in the 1000 realizations was relatively small, as shown by the error bars in Figure 6 denoting the standard deviation. For the setting shown in Figure 5 the mean computation time for VSTAP is 0.45 sec and for VARTA 1.55 sec with standard deviation 0.07 sec and 0.26 sec, respectively, indicating a significant difference even for small sample sizes.

Conclusion
The proposed method VSTAP can generate multivariate time series of arbitrary length with any given marginals and correlation structure, provided that the marginal distributions are continuous and the lagged correlation matrix is feasible. The general use of VSTAP lies in the linear piecewise approximation of the marginal transform from Gaussian to the given marginal, which allows for a closed form solution for the correlation transform from each component of the Gaussian lagged correlation matrix to the respective component of the given lagged correlation matrix. Thus any continuous marginal distribution, e.g. multi-modal or strongly skewed, can be sufficiently approximated at an accuracy depending on the number of breakpoints for the piecewise function. For all practical purposes, the accuracy converges with the number of breakpoints reaching the level of about 20, so that sufficient approximation of the sample marginal distribution is always obtained unless there are very few observed points, i.e. the time series is very short not allowing for the use of a sufficient number of breakpoints. Moreover, making use of the statistics of the joint doubly truncated Gaussian distribution, we could reach an analytic expression for the correlation transform. This allows for a straightforward and stable solution.
The VSTAP algorithm is also time effective, as the iterative scheme makes computation of a closed form expression for the correlation transform. This is to be compared to the numerical solution of the double integral form of the correlation transform used in the VARTA approach. We demonstrated with a simple example that piece-wise approximation results in much faster computations of the solution than numerical integration without any loss in accuracy. Still this may depend on the numerical integration scheme, which we did not investigate. So, besides the insight onto the correlation transform from Gaussian to target correlation provided by the closed-form approximation of the transform, a practical advantage of VSTAP is the derivation of the solution for the Gaussian correlation without the need of a time consuming two-dimensional numerical integration.
In VSTAP, we treated the problem of obtaining proper correlation matrices. Simple eigenvalue correction to make the correlation matrix positive semidefinite is not directly applicable to the lagged correlation matrix as it contains repeated entries. We introduced an iterative two-stage procedure that turned out to render positive definiteness in just few steps.
By construction VSTAP matches exactly the given marginals. In all simulations with different VAR processes, VSTAP could also match well the autoand cross-correlations for a sufficiently large number of lags when the length of the generated multivariate time series was up to moderately large, say up to about 4000. For larger lengths, some deviation could be observed for some specific auto-and cross-correlations, which however were rare and only in some of the studied systems (besides the presented simulation results for two systems, a number of other VAR systems of varying structure were tested). Thus for most practical purposes VSTAP generates proper multivariate time series that can be used for randomization tests or in stochastic simulation.
We have considered marginal transforms that deviate a lot from Gaussian and many different correlation structures. However, in our simulations, we have not encountered non-feasible correlation matrices as VSTAP would always provide sufficient solution, eventually after rendering positive definiteness of the lagged correlation matrix. It is therefore our intention to test VSTAP on special cases with more extreme marginals, e.g. positively skewed with a peak at zero, and stronger auto-and cross-correlations, e.g. typically expected from oscillating time series. The latter are often met in applications of the randomization test for nonlinearity, which is not discussed here but truly it has been the main motivation for this work. We leave this discussion and comparison to frequency-based methods, such as IAAFT for multivariate time series, including also time series from nonlinear dynamical systems, to future work. (h)