SIMULATING BIVARIATE STATIONARY PROCESSES WITH SCALE-SPECIFIC CHARACTERISTICS

By modifying and generalizing the wavelet-based approach of approximately simulating univariate long-memory processes that is available in the literature, we propose a methodology for simulating a bivariate stationary process, whose components exhibit different relationships at different scales. We derive the formulas for the autocovariance and cross-covariance sequences of the simulated bivariate process. We provide a setting for the parameters of the simulation which might generate a bivariate time series resembling that of stock log returns. Using this setting, we study the properties of our methodology via Monte Carlo simulation.


Introduction
In univariate analysis, wavelets are useful for studying the dynamics of a time series as a function of scale and time. In our paper, we will be predominantly interested in the scale aspect. A scale can be informally considered proportional to the reciprocal of frequency with long scales (i.e., low frequencies) being associated with the long-term behavior of the time series and short scales (i.e., high frequencies) with abrupt changes. A time series can generally exhibit different variability at different scales. For example, it can be highly variable at short scales with its dynamics consisting of prominent abrupt changes, jumps, etc. Alternatively, it can be highly variable at long scales exhibiting noticeable trends and long-term patterns. We will call the characteristics that capture the variability of a time series across scales scale-specifi c characteristics of variability. One well-known scale-specifi c characteristics of variability is wavelet variance (see, e.g., Percival and Walden, 2002, Chapter 8).
In case of bivariate analysis, wavelets are suitable for exploring scale-specifi c relationships, i.e., relationships that change across scales. Tools suitable for studying scale-specifi c relationships include wavelet covariance, wavelet correlation, wavelet cross-covariance and wavelet cross-correlation (see, e.g., Percival and Walden, 2002;Whitcher et al., 1999;Whitcher et al., 2000;Walden, 2000a and2000b). These tools have already been used in fi nance and economics. Ramsey and Lampart (1998a) and Ramsey and Lampart (1998b) studied the relationship between money and income and between expenditure and income, respectively. Atkins and Sun (2003) used wavelets to explore the relationship between interest rates and infl ation and Gençay et al. (2005) employed wavelets to explore relationships between the return of a portfolio and its beta. The results of these papers suggest that wavelets can be useful for revealing the nature of association between economic and fi nancial time series.
In this paper, we deal with a related problem, namely how to build models of univariate and bivariate time series with prescribed scale-specifi c characteristics of variability and scale-specifi c relationships. In the univariate case of generating a single time series with prescribed scale-specifi c characteristics of variability, we modify the approach of Percival and Walden (2002, Section 9.2), which builds upon Wornell (1993 and1996) and McCoy and Walden (1996). Afterwards, we suggest a way how this approach might be extended to the bivariate case of generating a pair of time series with prescribed scale-specifi c relationships. We fi nd the formulas for the traditional autocovariance and cross-covariance sequences of the simulated bivariate process.
Related methodologies for simulating time series with scale-specifi c characteristics have already appeared in the literature; see, for example, Nason et al. (2000) and Sanderson et al. (2010). The relation of these methodologies to our approach will be discussed in Section 7.
Our paper is organized as follows. In Section 2, we defi ne the discrete wavelet transform. In Section 3, we introduce the notion of wavelet variance, wavelet covariance and wavelet correlation. In Section 4, we modify the approach of Percival and Walden (2002, Section 9.2) and present a methodology for simulating univariate time series with scale-specifi c characteristics of variability. In Section 5, we generalize the methodology to the case of simulating bivariate time series with prescribed scale-specifi c relationships. In Section 6, Monte Carlo simulation is used to verify that our methodology works well. Section 7 compares our models to the models of Nason et al. (2000) and Sanderson et al. (2010). Section 8 concludes, pointing out the relevance of our approach to fi nancial time series.

The discrete wavelet transform
Our introduction to the discrete wavelet transform (DWT) is based on Percival and Walden (2002).

The DWT
Let us assume a real-valued time series {X t : t = 0, 1, …, N -1} of length N, where N = 2 J , J being a positive integer. This restriction on the length of the time series is necessary for defi ning the DWT. Let X be an N-dimensional vector consisting of the values of the time series {X t }, the fi rst element of X being X 0 , the second X 1 , etc. X can be considered a random vector since the time series can be considered a portion of length N of a real-valued discrete-time stochastic process.
The DWT of X is defi ned as where W is a real-valued N-dimensional vector of the DWT coeffi cients. Ω is an N × N non-random real-valued orthonormal matrix, which defi nes the DWT. The elements of Ω are fully determined if the length of the time series, N = 2 J , and the type of wavelets (e.g., Haar, Daubechies, coifl et etc.) used in the analysis are given. For example, for N = 2 3 = 8 and Haar wavelets, Ω is given as 1  1  1  1  0  0  0  0  2  2  2  2  1  1  1  1  0  0  0  0  2  2  2  2  1  1  1  1  1  1  1  1   8  8  8  8  8  8  8  8  1  1  1  1  1  1  1  1   8  8  8  8  8  8  8  8 Haar and the size of Θ J is 1 × N. The (r + 1)th (for r = 1, …, N j -1) row of Ω j is the rth row of Ω j circularly shifted by 2 j units. The partitioning of Ω implies the following partitioning of W into J + 1vectors W 1 , W 2 , …,W J , V J : where W j ≡ Ω j X (for j = 1, 2, …, J) is an N j -dimensional vector of the jth level wavelet coeffi cients, and V J = Θ J X is a one-dimensional vector consisting of a single Jth level scaling coeffi cient. The elements of W j (for j = 1, 2, …, J) are (for many types of wavelets) associated with changes of (weighted) averages of the input time series (X), these averages being calculated at scale 2 j-1 . V J , on the other hand, has a different meaning: it is associated with the average of the input time series calculated at scale 2 J .

The DWT and linear filtering
The DWT coeffi cients can also be obtained as a result of circular fi ltering the time series {X t : t = 0, 1, …, N -1} with a special set of J + 1 fi lters and via subsampling of the result subsequently. The J + 1 fi lters are denoted as {h 1,t : t = 0, ... , L 1 -1}, {h 2,t : t = 0, ... , L 2 -1}, ..., {h J,t : t = 0, ... , L J -1} and {g J,t : t = 0, ... , L J -1}, where L 1 , L 2 , ..., L J are the lengths of these fi lters. The length L j (for j = 2, ..., J) is related to the length L 1 as The exact values of these J + 1 fi lters are fully determined by the type of wavelets used in the analysis and are directly related to the elements of Ω.
Now, let us express W j (for j = 1, 2, …, J) and V J as where W j,t (for t = 0, …, N j -1) is the (t + 1)th element of W j and V J,0 is the only element of V J . We can write 1 , , 2 ( 1) 1 mod 0 , for 1,..., , 0,..., 1, where "mod" denotes the operation modulo. Equations 8 and 9 imply that the elements of W j (for j = 1, ..., J) and V J can be obtained via circular fi ltering of {X t } with {h j,t } and {g J,t }, respectively, and via the consequent subsampling of the result by 2 j and 2 J , respectively. Those elements of W j (for j = 1, ..., J) for which the omission of the modulo operation in Equation 8 would lead to the coeffi cients not being properly defi ned, are called boundary coeffi cients. The value of the boundary coeffi cients is affected by the assumption inherent in Equation 8, namely the fact that the input time series is treated as if it was periodic with the period equal to N, the values preceding X 0 being X N -1 , X N -2 etc. Since this assumption is often unrealistic in practical settings, boundary coeffi cients are often excluded from further analysis. As documented in Percival and Walden (2002, pages 136 and 146), the fi rst coeffi cients of W j (for j = 1, ..., J) are boundary coeffi cients (ceiling(b) is defi ned as the smallest integer greater than or equal to b). The remaining coeffi cients of W j are called non-boundary coeffi cients. Similarly, the coeffi cient V J,0 of Equation 9 will be treated as boundary in our paper if the modulo operation cannot be left out of its defi nition.

Inverse DWT
The DWT is an orthonormal transform which implies that the inverse DWT can be written as Making use of the partitioning of Ω and W described in Section 2.1, the matrix multiplication of Equation 13 can be rewritten as Equation 14 implies that the input time series (X) can be synthesized from the DWT coeffi cients (W). The term Ω j T W j might be considered the contribution of scale 2 j-1 to this synthesis.

Wavelet covariance, variance and correlation
In the next paragraphs, we introduce the notions of wavelet and scaling covariance (Section 3.1) and wavelet and scaling correlation (Section 3.4), which are helpful for exploring scale-specifi c relationships between a pair of processes. Wavelet and scaling variances that can be utilized for exploring scale-specifi c characteristics of variability of a single process are introduced in Section 3.3.

Wavelet and scaling covariance
Let us assume two real-valued time series {X t : t = 0, 1, …,N -1} and {Y t : t = 0, 1, …, N -1}, both of a length N = 2 J , J being a positive integer. Let these time series be portions of length N of two zero-mean components, {X t : t = … , -1, 0, 1, …} and {Y t : t = … , -1, 0, 1, …}, of a real-valued bivariate discrete-time stochastic process. Let the cross-covariance sequence between the components of the bivariate process be defi ned as where τ is the lag. Let the cross-spectrum between the components be defi ned as the Fourier transform of {s XY,τ : τ = … , -1, 0, 1, …}, i.e., where i is the imaginary unit and f is a real-valued non-random variable called frequency. In this paper, the cross-covariance sequence will always be assumed to be symmetric around zero. Under this assumption, the cross-spectrum is obviously a real--valued, even and periodic function of frequency with the period equal to unity. Taking the inverse Fourier transform of such a cross-spectrum, we get Setting τ = 0 in Equation 17, we obtain which implies that the cross-spectrum informs us about the contribution of different frequency bands to the covariance s XY,0 .

Cross-covariance matrix of the DWT coefficients
Let the values of the time series {X t : t = 0, 1, …, N -1} and {Y t : t = 0, 1, …, N -1} be stacked into two N-dimensional zero-mean vectors denoted as X and Y. The cross--covariance matrix Σ XY between X and Y is a symmetric Toeplitz matrix given as Let W X ≡ ΩX and W Y ≡ ΩY be the DWT coeffi cients of X and Y. The cross-covariance matrix between W X and W Y is given as Let us stack the diagonal of this matrix into an N-dimensional vector P, which can be partitioned into J + 1 vectors P 1 , P 2 , …, P J , Q J in the following way: where P j (for j = 1, …, J) is an N j -dimensional vector and Q J one-dimensional. Those elements of P j (for j = 1, …, J) which are associated with non-boundary coeffi cients of W X,j and W Y,j (where W X,j ≡ Ω j X and W Y,j ≡ Ω j Y) are obviously all the same and equal to C XY,j = 2 j η XY,j. If the only element of Q J is also associated with non-boundary coeffi cients, its value is obviously equal to B XY,J = 2 J ξ XY,J .
Following a similar line of argument to that of Percival and Walden (2002, p. 348), we can argue that two requirements need to be met so that the off-diagonal elements Ψ XY are approximately zero. Firstly, we have to require that {h 1,t }, {h 2,t }, ..., {h J,t } and {g J,t } be band-pass fi lters for different (i.e., non-overlapping) pass bands. This requirement is approximately fulfi lled as has already been stated in Section 2.2. Secondly, the cross-spectrum S XY (.) must be approximately constant in each of the following ranges of frequencies: Even though this assumption is not met in general, it might be assumed to be valid for many real-life bivariate processes.
The Jth level scaling correlation is defi ned as , , , ,

Simulation of univariate stationary processes with prescribed wavelet and scaling variances
We introduce a way of simulating univariate processes with prescribed wavelet and scaling variances. The core of the methodology is given in Section 4.1, being a modification of the approach of Percival and Walden (2002, Section 9.2), which is based on the work of Wornell (1993 and1996) and McCoy and Walden (1996). The modifi cation lies in the fact that while Percival and Walden (2002, Section 9.2) aim at approximately simulating a long-memory process with a given spectrum using wavelets, we aim at exactly simulating a stationary process with prescribed wavelet and scaling variances. As a part of this modifi cation, we have to solve a system of linear equations for a set of free parameters (see Section 4.3) to ensure that the simulated process has exactly the prescribed characteristics. We also derive a formula for the autocovariance sequence of the simulated process (Section 4.2), emphasizing that it consists of contributions of different scales.

Modification of the approach of Percival and Walden
Ψ XX is approximately diagonal for many types of processes (see Section 3.3). Even though the uncorrelatedness of the components of W does not generally imply their independence, it can sometimes be the case (for example if W has a multivariate normal distribution). This offers an appealing way of simulating stationary processes with prescribed wavelet and scaling variances. More specifi cally, let us simulate a portion, denoted as {X t : t = 0, …, 4N -1} (of a length 4N = 2 J+2 , J being a positive integer) of a zero-mean stationary process with wavelet variances η XX,1 , η XX,2 , …, η XX,J+2 and the (J + 2)th level scaling variance ξ XX,J+2 . Even though we simulate a portion of a length 4N in the fi rst place, we will eventually use only the fi rst N simulated values due to a reason which becomes obvious at the end of Section 4.4.
Step 2: The inverse DWT of Equation 13 can be applied to U, the result being Ω T U.
Step 3: Ω T U encodes a portion of length 4N of a stochastic process, which we, among others, require to be stationary. However, this requirement is not met, because the autocovariance sequence of the portion of the stochastic process encoded in Ω T U is not constant in time (see Percival and Walden, 2002, p. 356). Percival and Walden (2002, p. 356) argue that stationarity can be achieved by applying a random circular shift to Ω T U. More specifi cally, let the 4N × 4N matrix Π be defi ned in the following way: The result of the random circular shift applied to Ω T U is a 4N-dimensional vector defi ned as where k denotes a random variable distributed uniformly over integers 0, 1, …, 4N -1. It follows that the covariance matrix of X is given as where Σ UU is the covariance matrix of U and circ is such an operator whose input can be any (4N × 4N) matrix Δ, the output being defi ned as

Obtaining the exact values of wavelet and scaling variances
Apart from being stationary, we require the simulated process to have wavelet variances equal to η XX,1 , η XX,2 , …, η XX,J+2 and the (J + 2)th level scaling variance equal to ξ XX,J+2 . Making use of Equation 36 and the partitioning of Ω given in Equation 3, we can obtain 2 2 , , 2 which implies that Σ XX is a sum of several matrices each associated with a different scale. Further, let Ω i X (for i = 1, …, J + 2) be the vector of the ith level wavelet coeffi cients of X and let Γ i = Ω i Σ XX Ω i T be the covariance matrix of Ω i X. It follows that 2 2 , , 2 2 2 1 2 ( ) 2 ( ) . Equation 36) and the (r + 1)th (for r = 1, …, 4N/2 i -1) row of Ω i is the rth row of W i circularly shifted by an amount 2 i (see Section 2.1). Consequently, the diagonal of Γ i (for i = 1, …, J + 2) is constant and is determined by one single number, which we require to be equal to 2 i η XX,i . Now, there are (J + 3) parameters α XX,1 , α XX,2 , …, α XX,J+2 and β XX,J+2 , which determine the value of the diagonal of Γ i (see Equation 39). However, the total number of matrices Γ i is J + 2, because i runs over integers from 1 through J + 2. Moreover, there is one additional covariance matrix, corresponding to the (J + 2)th level scaling coeffi cients of X, whose diagonal is constant and required to be equal to 2 J+2 ξ XX,J+2 . All together, we thus have (J + 3) linear equations for (J + 3) unknown parameters α XX,1 , α XX,2 , …, α XX,J+2 and β XX,J+2 . If we solve for these parameters and if all are non-negative, then they are valid inputs into step 1 of Section 4.1 and lead to exact simulation with the required wavelet and scaling variances.

Symmetry of the autocovariance sequence
X might be considered a realization of a stationary stochastic process {X t : t = 0, …, 4N -1}. Since Σ XX is a circulant matrix, the autocovariance sequence of {X t : t = 0, …, 4N -1} is symmetric around the lag τ = 4N/2. To avoid this undesirable property, we apply the suggestion of Percival and Walden (2002, p. 357) and use only the fi rst N simulated values (out of the total number 4N) for the fi nal output. Thus, if the simulated time series is required to have the length N at the output, we have to simulate 4N values in the fi rst place.

Simulation of bivariate stationary processes with prescribed wavelet and scaling variances/correlations
In this section, we propose a generalization of the methodology presented in Section 4 to the bivariate case. The two components that constitute the bivariate process that we want to simulate will comove, the nature of the comovement being different on different scales. By construction, the cross-covariance sequence between the components will be symmetric around zero.

Methodology
Let us simulate portions {X t : t = 0, …, 4N -1} and {Y t : t = 0, …, 4N -1}, both of a length 4N = 2 J+2 , of two zero-mean components of a real-valued bivariate stationary process. Even though 4N values of each component are simulated in the fi rst place, we will use only the fi rst N simulated values in the end (because of the argument of Section 4.4). Let the jth level (for j = 1, …, J + 2) wavelet variance of the fi rst and second component (of the bivariate process) be denoted as η XX, j and η YY,j , respectively. Similarly, let the (J + 2)th level scaling variances of the components be denoted as ξ XX,J and ξ YY,J , respectively. Moreover, let the fi rst through the (J + 2)th level wavelet correlations and the (J + 2)th level scaling correlation between the components be denoted as ρ XY,1 , ρ XY,2 , …, ρ XY,J and φ XY,J .
Let the two portions of a length 4N of the components of the bivariate process be obtained by the approach of Section 4.1 and stacked into two 4N-dimensional vectors X and Y. We may write where the realization of the random variable k is the same in both the equations. U and V are two 4N-dimensional zero-mean generating vectors. U is formed analogously to the way presented in Section 4.1, i.e., its elements are 4N independent random variables with the unique variances 2 1 α XX,1 , 2 2 α XX,2 , …, 2 J+2 α XX,J+2 and 2 J+2 β XX,J+2 set up in such a way so that the jth level (for j = 1, …, J + 2) wavelet variances and the (J + 2)th level scaling variance of the fi rst component are equal to η XX,1 , η XX,2 , …, η XX,J+2 and ξ XX,J+2 . V is constructed analogously, with the unique variances 2 1 α YY,1 , 2 2 α YY,2 , …, 2 J+2 α YY,J+2 and 2 J+2 β YY,J+2 of its elements set up in such a way so that the jth level (for j = 1, …, J + 2) wavelet variances and the (J + 2)th level scaling variance of the second component are equal to η YY,1 , η YY,2 , …, η YY,J+2 and ξ YY,J+2 .
The second-order relationship between X and Y is described by the following cross-covariance matrix (see also Section 3.2): where   .
To ensure that the cross-covariance sequence between the components (of the bivariate process) is symmetric around zero, Σ UV has to be symmetric. Moreover, based on the argument at the end of Section 3.2, the off-diagonal elements of Σ UV can be set up equal to zero without losing much generality. As a result, we assume Σ UV to be diagonal.

Obtaining the exact values of wavelet and scaling correlation
Following a similar line of argument to that of Section 4.2, we must set up the (J + 3) values of the parameters γ XY,1 , γ XY,2 , …, γ XY,J+2 and δ XY,J+2 in such a way that the wavelet correlations between the components are equal to the required values ρ XY,1 , ρ XY,2 , …, ρ XY,J+2 and the (J + 2)th level scaling correlation to φ XY,J+2 . This can be achieved by solving an appropriate system of (J + 3) linear equations for the values of the parameters γ XY,1 , γ XY,2 , … , γ XY,J+2 and δ XY,J+2 . Since these parameters play the role of correlation coeffi cients, their values must lie in the range from -1 (included) to 1 (included) to be valid inputs into our methodology. Shall the solution to the system of linear equation not satisfy this constraint, it is not possible to simulate a bivariate process with the required wavelet and scaling variances and correlations.

Illustration
Using our methodology, we generate two time series that might resemble (in some sense) time series of stock log returns, which are strongly negatively correlated on short scales and positively on long scales.

Setting of the simulation
We simulate a portion of a bivariate stationary stochastic process with zero-mean components. We use Haar wavelets and choose J = 9, which implies that 4N = 2 J+2 = 2048. Only the fi rst N simulated values will be used for the fi nal output. The univariate distribution of individual elements of the generating vectors U and V is normal.
The components of the simulated bivariate process are required to have the jth level wavelet variances (for j = 1, …, 11) and the 11th level scaling variance in agreement with Table 1, which reports the 2 j th multiples of these characteristics. These scale--specifi c characteristics of variability are chosen in such a way so that both the components exhibit white noise behavior, the spectrum of both processes being constant and equal to unity. Consequently, the wavelet and scaling variances of both the processes are calculated based on the theory and formulas of Sections 3.3 and 3.1. Moreover, Table 1 also reports the jth level wavelet correlations (for j = 1, …, 11) and the 11th level scaling correlation between the components. These characteristics of scale--specifi c relationships indicate a negative relationship on short scales and a positive one on long scales. For our setting of the parameters, we can also write ρ XY,j = 2 j η XY,j (for j = 1, …, 11) and φ XY,11 = 2 11 ξ XY,11 .

Monte Carlo simulation
To verify that our approach simulates the bivariate process with the prescribed scale--specifi c characteristics correctly, we run the simulation with the setting of Section 6.1 5000 times. For each simulation wavelet and scaling variances of the fi rst component, wavelet and scaling covariances between the components, the autocovariance sequence of the fi rst component and the cross-covariance sequence between the components are estimated from the fi nal output of length N = 2 9 = 512. The reason for not estimating the univariate characteristics of the second component stems from the fact that these are identical to those of the fi rst component.
More specifi cally, the unbiased point estimates of 2 j η XX,j (for j = 1, …, 9), 2 9 ξ XX,9 , 2 j η XY,j (for j = 1, …, 9) and 2 9 ξ XY,9 are calculated for each simulation. Afterwards, these estimates are averaged across the 5000 simulations. These averages are subtracted from the prescribed (theoretical) values, 2 j η XX,j (for j = 1, …, 9), 2 9 ξ XX,9 , 2 j η XY,j (for j = 1, …, 9) and 2 9 ξ XY,9 , and the differences are reported in Table 3 together with the estimates of the standard errors of the averages. The reason for estimating wavelet and scaling covariances rather than wavelet and scaling correlations is the fact that we can easily construct an unbiased estimate of wavelet and scaling covariance, which is not the case with estimating wavelet and scaling correlation 1 .
The averages of estimates of 2 9 ξ XX,9 and 2 9 ξ XY,9 are reported in the last column of Table 3, being subtracted from the prescribed (theoretical) values of 2 9 ξ XX,9 and 2 9 ξ XY,9 . However, the prescribed values of 2 9 ξ XX,9 and 2 9 ξ XY,9 are not available at the initial set up of the simulation (see Table 1). Fortunately, we can obtain 2 9 ξ XX,9 and 2 9 ξ XY,9 in terms of η XX,10 ,η XX,11 ,ξ XX,11 ,η XY,10 ,η XY,11 and ξ XY,11 (of Table 1), which cannot be estimated directly since the fi nal output of the simulation is of a length N rather than 4N. More specifi cally, it follows from Equation 28 and Section 3.3 that Similarly, we calculate the sample autocovariance sequence of the fi rst component and the sample cross-covariance sequence 2 for each of the 5000 simulations. Then, we calculate the average of the sequences across the 5000 simulations. These averages are subtracted from the theoretical sequences of Figure 1, the differences at individual lags being plotted in Figure 2 together with the estimates of the standard errors of the averages. Again, the results suggest that the true autocovariance and the true cross--covariance sequence of the simulated bivariate process are indeed identical to the theoretical ones of Figure 1. Left-hand plot: The difference (thin line) between the average of the 5000 sample autocovariance sequences from the Monte Carlo simulation and the theoretical autocovariance sequence of Figure 1. Estimates of the standard errors of the average are plotted with the thick line. Right-hand plot: Analogously for the cross-covariance sequence.

Examining one realization of the bivariate process
The left-hand plot of Figure 3 presents one realization of the bivariate process of the fi nal length N = 2 9 . Point estimates of 2 j η XX,j (for j = 1, …, 9) and 2 9 ξ XX,9 and of 2 j η XY,j (for j = 1, …, 9) and 2 9 ξ XY,9 are presented in Table 4. The sample autocovariance sequences of both the components and the sample cross-covariance sequence are presented in Figure 4.
To appreciate the importance of the positive relationship between the components on long scales, it is instructive to plot the cumulative sum of the values of both the components -see the right-hand plot of Figure 3. The fact that the components (before they were cumulatively summed) appear to be simultaneously negatively correlated (as judged by the sample cross-covariance sequence of Figure 4) might suggest that the time series of the cumulative sums should move in opposite directions. However, this is true only for short time horizons, while in the long run these time series share the same direction, which is in agreement with the setting of the simulation. Table 4 Point estimates of scale-specifi c characteristics of the simulated bivariate process calculated for the realization in the left-hand plot of Figure 3. Left-hand plot: Simulation of the bivariate process with the prescribed scale-specifi c characteristics of Table 1. For the plotting purposes only, the realization of the second component is shifted by 7 units along the y-axis towards minus infi nity. Right-hand plot: Sequences of cumulative sums.

Figure 4
The sample autocovariance sequence of both components of the simulated bivariate process (top row, solid line) and the sample cross-covariance sequence between the components (bottom row, solid line). Thin vertical lines are plotted at the lag 0 and thin horizontal lines at the value of the autocovariance/cross-covariance equal to 0.

Comparison of our model with the literature
Equation 13 represents X in terms of W. Since Ω is orthogonal and thus invertible (Section 2.3), this representation is unique in the sense that, if X is given, the vector W, which satisfi es Equation 13, is unique and given by Equation 1. Elements of W can be considered coordinates of X with respect to an orthonormal basis, the basis vectors being the columns of Ω T and thus being associated with particular scales. The approach of the univariate simulation of Section 4 thus consists in reconstructing X from its coordinates (ΩΠ k Ω T U) with respect to the orthonormal basis (see Equation 35).
Since the main diagonal of Ω j T Ω j is generally not constant, the univariate distribution of a given element of X is that of a mixture of Gaussians.
There are other fl avors of wavelet transform such as the MODWT (see Percival and Walden, 2002), which projects X onto an overcomplete basis (frame), whose basis vectors are again associated with particular scales. While the use of the overcomplete basis leads to a non-unique representation of X in terms of the coordinates of X with respect to the overcomplete basis, its advantage lies in the fact that it provides a better way (when compared to the orthonormal basis of the DWT) of exploring "localized" structures in X, which can be useful for studying non-stationary time series. This is the reasoning of Nason et al. (2000), who proposed a univariate locally stationary wavelet model based on an overcomplete basis.
The model of Nason et al. (2000) is primarily a model for univariate non-stationary time series with a time-varying autocovariance sequence. It can also be used to model univariate stationary time series with scale-specifi c characteristics of variability and for the simulation thereof. In that case, the simulated time series is reconstructed from the coordinates in the overcomplete basis; see Equation 7 of Nason et al., 2000. It follows from this equation that the simulated time series is necessarily Gaussian (as opposed to the non-Gaussian case of our model) if the joint distribution of the coordinates is multivariate normal. The autocovariance sequence of the simulated time series can be decomposed into so-called autocorrelation wavelets (see Nason et al., 2000, Section 2.3 and Equation 14). This decomposition is analogous to our decomposition of Equations 38 and 40.
Our simulation of the bivariate time series of Section 5 is a generalization of the univariate simulation of Section 4. Similarly, the model of Sanderson et al. (2010) is an extension of the locally stationary wavelet time series model of Nason et al. (2000) to the bivariate case. It aims at modeling non-stationary time series with a time-varying autocovariance and cross-covariance structure.
The models of Nason et al. (2000) and Sanderson et al. (2010) are supported with a well-developed theory and are generally more fl exible when compared to our models of Section 4 and Section 5, since they can better handle time-varying second-order structures. On the other hand, since our models start from an orthonormal basis rather than an overcomplete one, their construction is more intuitive and may be more readily understood. The non-Gaussian nature of the simulated time series may also be convenient in some situations. Moreover, as a part of our models, we have also provided a way of obtaining simulated time series with exact scale-specifi c characteristics of variability and relationships.

Conclusion and discussion
We have modifi ed and generalized the wavelet-based approach of Percival and Walden (2002), which was originally used for approximate simulation of univariate long--memory processes. More specifi cally, we proposed a methodology for simulating bivariate stationary processes exhibiting scale-specifi c relationships. Our simulation is exact in the sense that it achieves the required scale-specifi c characteristics. This was accomplished by solving a system of linear equations for several free parameters of the simulation. Using Monte Carlo simulation, we have verifi ed that our methodology works correctly.
Our methodology might be suitable wherever scale-specifi c relationships are assumed to be present, which might possibly not be a scarce situation in fi nance and economics. This belief is supported by the fact that fi nancial markets are formed by several groups of traders (market makers, intraday traders, hedge funds, portfolio managers, investment banks, etc.), which employ different tools of market analysis and have different information at hand, different risk profi les and trading preferences and constraints. It can be argued that such differences among the groups might imply different time scales on which each group operates. Since each time scale is thus predominantly "operated" by a different group of traders with different characteristics, it can be expected that such a fi nancial market, called the heterogeneous market, might exhibit different behavior on different scales (see also Zumbach and Lynch, 2001;Müller et al., 1993;Müller et al., 1997).