Matching Kollo measures

Abstract

We motivate the advantages of using the Kollo measures, relative to other types of third and fourth moments of multivariate systems, and explore their Monte Carlo simulation and bootstrapping errors. Then we derive necessary and sufficient conditions for simultaneously matching any given mean vector, covariance matrix, Kollo skewness, and Kollo kurtosis. The specification of a suitable orthonormal basis greatly simplifies these moment-matching conditions. We offer semi-closed-form solutions to increase computational efficiency. In this respect, we compare our approach to two competing methods, which anyway can only match Kollo skewness and not the kurtosis at the same time. Ours is the first method for exactly matching all four multivariate moments simultaneously.

Keywords:

• Simulation
• statistics
• finance

1. Introduction

Consider an n-dimensional random variable ${\mathbf{X}}_n=\left[X_1,\dots ,X_n\right]\prime$ with a mean vector ${\mathit{\mu }}_n=\mathbb{E}[{\mathbf{X}}_n]$ and a covariance matrix ${\mathbf{\Sigma }}_n=\mathbb{V}[{\mathbf{X}}_n].$ Mitchell (1989) introduced the third and fourth moments as the skewness tensor $\mathbb{E}[{\mathbf{Y}}_n\otimes {\mathbf{Y}}_n^{\prime }\otimes {\mathbf{Y}}_n]$ and kurtosis tensor $\mathbb{E}[{\mathbf{Y}}_n\otimes {\mathbf{Y}}_n^{\prime }\otimes {\mathbf{Y}}_n\otimes {\mathbf{Y}}_n^{\prime }]$ where ${\mathbf{Y}}_n$ is the standardized version of ${\mathbf{X}}_n,$ i.e., ${\mathbf{Y}}_n={\mathbf{\Sigma }}_n^{-1/2}({\mathbf{X}}_n-{\mathit{\mu }}_n)$ and ⊗ denotes the Kronecker product. In other words, the elements in the skewness tensor take the form $\mathbb{E}[Y_i{Y}_j{Y}_k]$ and those of the kurtosis tensor take the form $\mathbb{E}[Y_i{Y}_j{Y}_k{Y}_t],$ for $i,j,k,t=1,\dots ,n.$ The tensors contain the third and fourth order central moments and all the co-moments of the marginal variables, which include $n(n+1)(n+2)/6$ and $n(n+1)(n+2)(n+3)/24$ distinct elements, respectively. The number of elements increases rapidly with n. For that reason, these full tensor moments are impractical for use in large dimensions and several more parsimonious representations of multivariate skewness and kurtosis have been proposed.

Many authors advocate the need for models to match the information contained in measures of skewness and kurtosis, for instance in finance (Ryoo, 2007), biotech (Kuruvilla & Gunavathi, 2014), and engineering (Love et al., 2013). Indeed, moment-matching techniques have been applied in a variety of disciplines: Heath et al. (2018) use moment-matching simulation for probabilistic sensitivity analysis in economics; Pellegrino and Sabino (2014) apply moment-matching to reduce the dimension of the pricing problem in finance; and Ching et al. (2009) improve the accuracy of the estimation by using moment-matching in engineering. Yet none can help match even the mean vector and covariance matrix exactly, let alone the higher moments.

Some moment-matching techniques are specifically designed to reduce the sampling error in multivariate moments. In this strand of the literature both Bardsley (2017), Høyland et al. (2003), and Looney (1983) propose a method to generate data based on moments; Foldnes and Olsson (2016) introduce a simple copula structure to target the mean vector, the covariance matrix and both the skewness and the kurtosis of each marginal; and Qu et al. (2020) propose a method to target the Mardia (1970) skewness and kurtosis. However, none of these use the tensor measures defined above, because they are intractable for practical applications except when n is small.

Many different summary measures have been advocated to reduce the size of the skewness and kurtosis tensors, including those proposed by Koizumi et al. (2014), Kollo (2008), Mardia (1970), and Móri et al. (1994). Several other measures are not even directly related to the full tensors.¹ All these measures are listed in chronological order in Table 1 and each column refers to a property that is not fulfilled by the measure.

Table 1. Weaknesses of popular skewness and kurtosis measures for multivariate variables.

Important properties	Tensor- based	All cross moments included	Non-scalar measure	Affine invariant
Mardia (1970)		✗^†	✗
Skewness and Kurtosis
Malkovich and Afifi (1973)	✗	Not applicable	✗
Skewness and Kurtosis
Isogai (1982)^‡	✗	Not applicable	✗	✗
Skewness
Srivastava (1984)	✗	Not applicable	✗	✗
Skewness and Kurtosis
Móri et al. (1994)		✗
Skewness and Kurtosis
Song (2001)	✗	Not applicable	✗
Kurtosis
Balakrishnan et al. (2007)	✗	Not applicable
Skewness and Kurtosis
Kollo (2008)
Skewness and Kurtosis
Koizumi et al. (2014)		✗^†	✗
Skewness and Kurtosis

^† All cross-moments are considered for the skewness but not kurtosis.

^‡ Applicable for distributions without finite first fourth moments.

First note that most of the measures that are not based on the full tensors require computation of the multivariate probability density function, which is rather complex, and cannot be used to establish solvable systems of moment-matching constraints for simulations targeting skewness and kurtosis.² Secondly, of the four methods that are derived from the full tensor moments, all but three ignore information by excluding some of the cross moments. This includes the popular Mardia (1970) skewness and kurtosis, which are easy to use but—like several other measures—are only scalar quantities. As such, they are unable to carry as much information about the full multivariate density as vector or matrix measures.

From Table 1 it is clear that the only measures that suffer from none of these drawbacks are those introduced by Kollo (2008). In the form of an $n\times 1$ vector and an n × n symmetric matrix, respectively, they may be written: (1a) $$\mathit{\tau }\left({\mathit{X}}_n\right)=\left[\begin{array}{c}\sum _{i,k}^n\mathbb{E}\left[Y_i{Y}_1{Y}_k\right],\dots ,\sum _{i,k}^n\mathbb{E}\left[Y_i{Y}_n{Y}_k\right]\end{array}\right]\prime ,$$(1a) (1b) $$\mathit{K}\left({\mathit{X}}_n\right)=\left[\begin{array}{ccc}\sum _{i,k}^n\mathbb{E}\left[Y_i{Y}_1{Y}_1{Y}_k\right]& \cdots & \sum _{i,k}^n\mathbb{E}\left[Y_i{Y}_1{Y}_n{Y}_k\right]\\ ⋮& \ddots & ⋮\\ \sum _{i,k}^n\mathbb{E}\left[Y_i{Y}_n{Y}_1{Y}_k\right]& \cdots & \sum _{i,k}^n\mathbb{E}\left[Y_i{Y}_n{Y}_n{Y}_k\right]\end{array}\right].$$(1b)

Unlike tensor moments, the Kollo measures are easy to compute for applications. But so are other measures, such as those of Mardia and Móri et al.³ In Appendix A, we provide some real-world examples that motivate the use of Kollo skewness and kurtosis in a practical setting, rather than the measures introduced by Mardia. Kollo measures also have many theoretical advantages. By retaining so many co-moments the Kolo measures are the most tractable attempt to retain as much information as possible from the full skewness and kurtosis tensors.⁴ The Mardia and Koizumi et al. measures suffer from a substantial loss of information because they are scalars. For example, Mardia skewness can only be positive, which cannot show the data is negatively skewed or positively skewed. A similar issue of information loss is shared by the Móri et al. measures which are a multivariate extension of Mardia measures but still omit many third and fourth order co-moments. Moreover, like the tensors themselves and several other related measures, the Kollo measures are invariant under affine transformations. Later we shall rely on this property to generate repeated simulations that exactly match some target values for Kollo skewness and kurtosis. This invariance property also enables concatenations of repeated samples. That is, the concatenated samples have the same moments as the sub-samples.

Our exact moment-matching methods are achieved by creating an orthonormal basis, the span of which includes a particular subspace containing all samples that have the exact Kollo skewness vector and Kollo kurtosis matrix. The benefit of using such a set of orthonormal vectors is huge: mathematically, matching the sample moments exactly means solving a large system of non-linear equations in which the number of unknowns greatly exceeds the number of constraints. The orthonormal basis also underpins a significant reduction in the size of the system to be solved, and we are able to derive analytic solutions that greatly enhance the computational efficiency of our method. It is also, to our knowledge, the only method for exactly matching both Kollo skewness and kurtosis in every simulation of the distribution.

Common techniques, such as Monte Carlo (MC) simulation (Hertz, 1964) or the statistical bootstrap (Efron, 1979) fail to ensure that each sample has the same moments as the target distribution. There is a large literature on their simulation errors in covariance matrices—see Jobson and Korkie (1980), Ledoit and Wolf (2004), and many others—but no previous study on MC or bootstrap simulation errors in the Kollo measures. Section 2 motivates the need to target higher moments exactly by demonstrating that the Kollo moment sampling errors from standard Monte Carlo and the statistical bootstrap can be very significant. This well-known property of sampling error is unlikely to be resolved completely by increasing the number of simulations alone (Ferraz & Moura, 2012; Oberkampf et al., 2002) and variance-reduction techniques (Saliby & Paul, 2009) come at the expense of model complexity.

Then Section 3 introduces the Ledermann et al. (2011) Random Orthogonal Matrix (ROM) simulation method, which encompasses MC and the bootstrap (as special cases of the L matrices which underpin this approach) and was originally developed as a simulation method for matching exact covariance matrices. Two different extensions of the ROM algorithm by Hanke et al. (2017) and Alexander et al. (2022) can target exact Kollo skewness by generating arbitrary values to pad the sample so that Kollo skewness matching conditions are solvable numerically. But a trial-and-error process is needed to find suitable arbitrary values and, even when such values are found, there is no analytical solution and its numerical solution cannot be guaranteed, so the algorithms can fail. Moreover, as Alexander et al. (2022) demonstrate, these methods risk generating simulated samples with undesirable marginal densities, depending on the arbitrary values that are found.

Having derived the necessary and sufficient conditions for a sample generated to match a target Kollo skewness and Kollo kurtosis in Section 3, Section 4 uses these conditions to derive theoretical results which first define an orthonormal basis, and then find analytic solutions for the coefficients which generate the required subspace of the span of this orthonormal basis. Here we also prove the invariance under sample concatenation. We conclude the paper with a discussion of our results in Section 5, in the context of the previous literature.

Our algorithms are written in Python and available on Github https://github.com/vivid1005/GKROM-algorithm.

2. Sampling errors in Kollo skewness and kurtosis

This section presents two experiments that demonstrate the need to use a simulation method that can target Kollo moments exactly. Let ${\mu }_i,{\sigma }_{\mathit{\text{ij}}},{\tau }_i,$ and κ_ij denote the elements of the target mean ${\mathit{\mu }}_n,$ covariance ${\mathbf{\Sigma }}_n,$ Kollo skewness ${\mathit{\tau }}_n,$ and Kollo kurtosis ${\mathbf{K}}_n,$ respectively, and let and ${\widehat{\mu }}_i,{\widehat{\sigma }}_{\mathit{\text{ij}}},{\widehat{\tau }}_i,$ and ${\widehat{\kappa }}_{\mathit{\text{ij}}}$ denote the corresponding elements of the simulated sample. Then the RMSEs for each moment are given by: (2) $$\sqrt{\sum _i^n{({\mu }_i-{\widehat{\mu }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\sigma }_{\mathit{\text{ij}}}-{\widehat{\sigma }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},\sqrt{\sum _i^n{({\tau }_i-{\widehat{\tau }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\kappa }_{\mathit{\text{ij}}}-{\widehat{\kappa }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},$$(2) for the mean vector, covariance matrix, Kollo skewness vector, and Kollo kurtosis matrix, respectively. Now, for a given dimension n of the system of random variables and a given size m for the simulations, we generate a sample and compute the mean vector, covariance matrix, Kollo skewness vector, and Kollo kurtosis matrix and record the RMSEs using (2)(2) $$\sqrt{\sum _i^n{({\mu }_i-{\widehat{\mu }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\sigma }_{\mathit{\text{ij}}}-{\widehat{\sigma }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},\sqrt{\sum _i^n{({\tau }_i-{\widehat{\tau }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\kappa }_{\mathit{\text{ij}}}-{\widehat{\kappa }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},$$(2) . We repeat this many times and then record the average RMSE and its standard error.

Tables 2 and 3 tabulate the results of the first, bootstrapping experiment using two real data sets, viz. the U.S. stock sector returns data described above and a temperature data sample from the “Global and City Yearly Average Temperatures 1820–2015” data set available on Data.world. It includes the yearly average temperatures between 1849 and 2013 of the 11 cities in the United States (viz. Austin, Boston, Chicago, Kansas City, Las Vegas, Los Angeles, New Orleans, New York, Philadelphia, San Francisco, and Washington). In both cases n = 11 with the target mean vector, covariance matrix, Kollo skewness vector, and Kollo kurtosis matrix derived from (1) the first 1000 observations in the sample, from 1 December 2018 to 31 December 2021 for sector returns and (2)(2) $$\sqrt{\sum _i^n{({\mu }_i-{\widehat{\mu }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\sigma }_{\mathit{\text{ij}}}-{\widehat{\sigma }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},\sqrt{\sum _i^n{({\tau }_i-{\widehat{\tau }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\kappa }_{\mathit{\text{ij}}}-{\widehat{\kappa }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},$$(2) 166 observations from 1849 to 2013 for temperature data, respectively. The moments are reported in the upper part of the tables. Then we apply the statistical bootstrap to these samples, to generate random samples of size m, for $m=100,500,$ and 1000. For any given m, we generate 100,000 samples of this size. For each sample, we record the RMSEs calculated using (2)(2) $$\sqrt{\sum _i^n{({\mu }_i-{\widehat{\mu }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\sigma }_{\mathit{\text{ij}}}-{\widehat{\sigma }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},\sqrt{\sum _i^n{({\tau }_i-{\widehat{\tau }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\kappa }_{\mathit{\text{ij}}}-{\widehat{\kappa }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},$$(2) and then compute the average and standard error of these RMSEs over all 100,000 repetitions and report the results in the lower part of the tables.

Table 2. Sampling errors in multivariate moments for statistical bootstrapping financial data.

${\mathit{\mu }}_{11}\times {10}^4$	${\mathbf{\Sigma }}_{11}\times {10}^4$
$$\left[\begin{array}{c}-0.20\\ 4.82\\ 4.09\\ 7.72\\ 3.81\\ 5.79\\ 4.56\\ 11.19\\ 6.15\\ 4.54\\ 6.28\end{array}\right]$$	$$\left[\begin{array}{c}\begin{array}{ccccccccccc}5.39& 2.67& 2.79& 1.93& 1.28& 1.62& 3.23& 2.20& 1.76& 1.43& 1.90\\ 2.67& 2.38& 2.17& 1.68& 1.21& 1.45& 2.33& 1.97& 1.48& 1.35& 1.63\\ 2.79& 2.17& 2.41& 1.73& 1.23& 1.49& 2.46& 2.01& 1.53& 1.38& 1.70\\ 1.93& 1.68& 1.73& 2.05& 1.04& 1.33& 1.81& 2.18& 1.63& 1.07& 1.43\\ 1.28& 1.21& 1.23& 1.04& 1.24& 1.08& 1.32& 1.29& 1.07& 1.25& 1.28\end{array}\\ \begin{array}{ccccccccccc}1.62& 1.45& 1.49& 1.33& 1.08& 1.56& 1.59& 1.70& 1.29& 1.20& 1.35\\ 3.23& 2.33& 2.46& 1.81& 1.32& 1.59& 3.07& 2.12& 1.66& 1.50& 1.85\\ 2.20& 1.97& 2.01& 2.18& 1.29& 1.70& 2.12& 2.91& 1.97& 1.30& 1.68\\ 1.76& 1.48& 1.53& 1.63& 1.07& 1.29& 1.66& 1.97& 2.01& 1.08& 1.33\\ 1.43& 1.35& 1.38& 1.07& 1.25& 1.20& 1.50& 1.30& 1.08& 2.11& 1.77\\ 1.90& 1.63& 1.70& 1.43& 1.28& 1.35& 1.85& 1.68& 1.33& 1.77& 2.27\end{array}\end{array}\right]$$
${\mathit{\tau }}_{11}$	${\mathbf{K}}_{11}$
$$\left[\begin{array}{c}\begin{array}{c}-3.19\\ -1.83\end{array}\\ -1.13\\ -3.48\\ 4.61\\ -0.46\\ 3.02\\ 4.26\\ 0.15\\ -4.78\\ -6.05\end{array}\right]$$	$$\left[\begin{array}{ccccccccccc}83.74& 48.97& 22.19& 30.72& 25.95& 17.63& 37.65& 18.12& 3.05& 30.91& 43.42\\ 48.97& 61.41& 16.25& 23.32& 29.09& 13.40& 28.86& 12.56& 0.75& 35.04& 28.77\\ 22.19& 16.25& 26.53& 9.83& 16.79& 6.09& 14.42& 11.92& 5.96& 5.81& 12.47\\ 30.72& 23.32& 9.83& 35.57& 6.62& 11.31& 14.40& 4.47& 4.58& 17.88& 28.18\\ 25.95& 29.09& 16.79& 6.62& 61.87& 4.67& 33.53& 28.75& 6.26& -2.93& 3.36\\ 17.63& 13.40& 6.09& 11.31& 4.67& 25.55& 7.33& 1.97& 3.07& 14.98& 17.75\\ 37.65& 28.86& 14.42& 14.40& 33.53& 7.33& 58.60& 33.24& 7.93& -0.66& 28.78\\ 18.12& 12.56& 11.92& 4.47& 28.75& 1.97& 33.24& 39.38& 7.33& -12.28& 9.94\\ 3.05& 0.75& 5.96& 4.58& 6.26& 3.07& 7.93& 7.33& 17.80& -4.64& 5.58\\ 30.91& 35.04& 5.81& 17.88& -2.93& 14.98& -0.66& -12.28& -4.64& 58.90& 21.77\\ 43.42& 28.77& 12.47& 28.18& 3.36& 17.75& 28.78& 9.94& 5.58& 21.77& 65.31\end{array}\right]$$
m	RMSE				Ratio
	$\mathit{\mu }\times {10}^4$	$\mathbf{\Sigma }\times {10}^4$	$\mathit{\tau }$	K	$\mathit{\mu }$	$\mathbf{\Sigma }$	$\mathit{\tau }$	K
100	$\underset{(24.470)}{50.324}$	$\underset{(5.150)}{7.868}$	$\underset{(7.562)}{19.906}$	$\underset{(25.260)}{179.332}$	$\underset{}{2.057}$	$\underset{}{1.528}$	$\underset{}{2.632}$	$\underset{}{7.099}$
500	$\underset{(12.590)}{27.669}$	$\underset{(2.319)}{3.612}$	$\underset{(5.399)}{17.295}$	$\underset{(35.914)}{139.581}$	$\underset{}{2.198}$	$\underset{}{1.558}$	$\underset{}{3.203}$	$\underset{}{3.887}$
1000	$\underset{(9.916)}{23.828}$	$\underset{(1.568)}{2.535}$	$\underset{(5.281)}{12.758}$	$\underset{(42.166)}{104.109}$	$\underset{}{2.403}$	$\underset{}{1.617}$	$\underset{}{2.416}$	$\underset{}{2.469}$

The upper table reports the mean (${\mathit{\mu }}_{11}$), variance (${\mathbf{\Sigma }}_{11}$), Kollo skewness (${\mathit{\tau }}_{11}$), and kurtosis (${\mathbf{K}}_{11}$) defined in (1) derived from a sample of 1000 returns on each of n = 11 U.S. stock sector index. In the lower table, the first column states the sample size of each bootstrap from which we calculate the RMSEs for the mean, covariance, Kollo skewness, and kurtosis, using (2)(2) $$\sqrt{\sum _i^n{({\mu }_i-{\widehat{\mu }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\sigma }_{\mathit{\text{ij}}}-{\widehat{\sigma }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},\sqrt{\sum _i^n{({\tau }_i-{\widehat{\tau }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\kappa }_{\mathit{\text{ij}}}-{\widehat{\kappa }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},$$(2) . Repeating the sample and again computing the RMSE 10,000 times, we report here the average and standard error (in parentheses, below) for each moment in columns 2–5. The last four columns are the average RMSE divided by its standard error.

Table 3. Sampling errors in multivariate moments for statistical bootstrapping temperature data.

${\mathit{\mu }}_{11}\times {10}^4$	${\mathbf{\Sigma }}_{11}\hspace{1em}\times {10}^4$
$$\left[\begin{array}{c}20.04\\ 7.61\\ 10.24\\ 11.98\\ 17.45\\ 15.88\\ 20.36\\ 9.77\\ 12.08\\ 14.45\\ 12.13\end{array}\right]$$	$$\left[\begin{array}{c}\begin{array}{ccccccccccc}0.38& 0.26& 0.28& 0.32& 0.12& 0.08& 0.29& 0.24& 0.23& 0.03& 0.23\\ 0.26& 0.73& 0.63& 0.54& 0.10& 0.12& 0.28& 0.67& 0.59& 0.11& 0.56\\ 0.28& 0.63& 0.81& 0.71& 0.07& 0.06& 0.31& 0.63& 0.57& 0.06& 0.58\\ 0.32& 0.54& 0.71& 0.76& 0.13& 0.10& 0.32& 0.53& 0.49& 0.08& 0.50\\ 0.12& 0.10& 0.07& 0.13& 0.37& 0.31& 0.06& 0.08& 0.06& 0.23& 0.06\\ 0.08& 0.12& 0.06& 0.10& 0.31& 0.32& 0.04& 0.10& 0.08& 0.25& 0.07\\ 0.29& 0.28& 0.31& 0.32& 0.06& 0.04& 0.33& 0.28& 0.27& 0.01& 0.27\end{array}\\ \begin{array}{ccccccccccc}0.24& 0.67& 0.63& 0.53& 0.08& 0.10& 0.28& 0.64& 0.58& 0.09& 0.56\\ 0.23& 0.59& 0.57& 0.49& 0.06& 0.08& 0.27& 0.58& 0.54& 0.07& 0.52\\ 0.03& 0.11& 0.06& 0.08& 0.23& 0.25& 0.01& 0.09& 0.07& 0.25& 0.06\\ 0.23& 0.56& 0.58& 0.50& 0.06& 0.07& 0.27& 0.56& 0.52& 0.06& 0.52\end{array}\end{array}\right]$$
${\mathit{\tau }}_{11}$	${\mathbf{K}}_{11}$
$$\left[\begin{array}{c}1.86\\ 1.19\\ -0.84\\ 2.26\\ 2.94\\ 0.42\\ 0.09\\ 1.48\\ -1.10\\ -0.09\\ 0.90\end{array}\right]$$	$$\left[\begin{array}{c}\begin{array}{ccccccccccc}17.45& 4.80& -3.33& 3.65& 7.94& -2.04& 1.88& 3.45& -2.60& 2.15& 7.39\\ 4.80& 14.27& -1.42& 3.30& 4.32& -0.65& 4.16& 4.08& -0.17& 1.89& 6.24\\ -3.33& -1.42& 14.79& -2.28& -2.75& 4.30& 0.95& -1.84& 3.78& -0.21& -3.14\\ 3.65& 3.30& -2.28& 12.26& 4.51& -0.18& 1.51& 2.22& -0.70& 1.97& 4.01\\ 7.94& 4.32& -2.75& 4.51& 17.81& -2.75& 4.79& 3.60& -2.83& 0.47& 6.56\end{array}\\ \begin{array}{ccccccccccc}-2.04& -0.65& 4.30& -0.18& -2.75& 14.20& 0.72& -0.36& 3.51& 1.78& -3.02\\ 1.88& 4.16& 0.95& 1.51& 4.79& 0.72& 12.51& 1.50& 2.23& 1.33& 2.23\\ 3.45& 4.08& -1.84& 2.22& 3.60& -0.36& 1.50& 11.79& -0.44& 1.28& 4.96\\ -2.60& -0.17& 3.78& -0.70& -2.83& 3.51& 2.23& -0.44& 11.17& 1.13& -2.24\\ 2.15& 1.89& -0.21& 1.97& 0.47& 1.78& 1.33& 1.28& 1.13& 11.79& 1.23\\ 7.39& 6.24& -3.14& 4.01& 6.56& -3.02& 2.23& 4.96& -2.24& 1.23& 17.38\end{array}\end{array}\right]$$
m	RMSE				Ratio
	$\mathit{\mu }\times {10}^4$	$\mathbf{\Sigma }\times {10}^4$	$\mathit{\tau }$	K	$\mathit{\mu }$	$\mathbf{\Sigma }$	$\mathit{\tau }$	K
100	$\underset{(0.100)}{0.215}$	$\underset{(0.342)}{0.663}$	$\underset{(1.333)}{5.141}$	$\underset{(8.617)}{23.280}$	$\underset{}{2.160}$	$\underset{}{1.939}$	$\underset{}{3.856}$	$\underset{}{2.702}$
500	$\underset{(0.045)}{0.096}$	$\underset{(0.150)}{0.293}$	$\underset{(0.705)}{2.147}$	$\underset{(5.740)}{8.771}$	$\underset{}{2.149}$	$\underset{}{1.951}$	$\underset{}{3.047}$	$\underset{}{1.528}$
1000	$\underset{(0.032)}{0.068}$	$\underset{(0.107)}{0.209}$	$\underset{(0.452)}{1.485}$	$\underset{(3.236)}{5.766}$	$\underset{}{2.126}$	$\underset{}{1.963}$	$\underset{}{3.288}$	$\underset{}{1.782}$

The upper table reports the mean (${\mathit{\mu }}_{11}$), variance (${\mathbf{\Sigma }}_{11}$), Kollo skewness $({\mathit{\tau }}_{11}),$ and kurtosis $({\mathbf{K}}_{11})$ defined in (1) derived using a sample of 1000 average yearly temperatures of n = 11 U.S. cities. In the lower table, the first column states the sample size of each bootstrap from which we calculate the RMSEs for the mean, covariance, Kollo skewness, and kurtosis, using (2)(2) $$\sqrt{\sum _i^n{({\mu }_i-{\widehat{\mu }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\sigma }_{\mathit{\text{ij}}}-{\widehat{\sigma }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},\sqrt{\sum _i^n{({\tau }_i-{\widehat{\tau }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\kappa }_{\mathit{\text{ij}}}-{\widehat{\kappa }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},$$(2) . Repeating the sample and again computing the RMSE 10,000 times, we report here the average and standard error (in parentheses, below) for each moment in columns 2–5. The last four columns are the average RMSE over its standard error.

As expected both the average RMSE and its standard deviation decrease with m. However, their ratio (shown in the right-hand panel) does not. The errors in moments that arise when using the bootstrap are particularly apparent when only 100 observations are randomly selected. When m = 100 the ratio of the average RMSE on the Kollo skewness to its standard deviation is 2.632 for the stock returns data and 3.856 for the temperature data, and for the Kollo kurtosis, the corresponding ratio is 7.099 and 2.702, respectively. But even when m = 1000, so we bootstrap 100,000 samples having the same size as the original one, the errors relative to the target moments derived from the original sample remain substantial (with ratios are 2.416 and 2.469 for the stock returns data and 3.288 and 1.782 for the temperature data on the Kollo skewness and kurtosis, respectively). It is worth noting that much of the research on sampling errors in simulation focuses on errors in the covariance matrix, but looking at the ratios in the right-hand panel it appears that the covariance matrix is the least inaccurate of the four moments.

Table 4 reports the results of the second experiment which assumes an n-dimensional $\mathcal{N}({\mathbf{0}}_n,{\mathbf{I}}_n)$ distribution for n = 5 or 10 and applies standard Monte Carlo simulation to generate samples of size m for $m={10}^3$ or 10⁶. As expected both the average RMSE and its standard error decrease as m increases. However, again their ratio (shown in the right panel) does not. It is clear that all simulations yield inaccurate values for the mean vector, covariance matrix, Kollo skewness vector, and Kollo kurtosis matrix. It also becomes more difficult to simulate accurate values as the dimension of the system increases, especially for the covariance matrix. This observation was used to motivate the introduction of moment matching methods (such as Ledermann et al., 2011) and Table 4 demonstrates that similar inaccuracies occur when using Monte Carlo simulation in large systems for which accurate third and fourth moments are also important.

Table 4. RMSE of multivariate moments in Monte Carlo simulation.

n	m	RMSE				Ratio
		$\mathit{\mu }\times {10}^2$	$\mathbf{\Sigma }\times {10}^2$	$\mathit{\tau }\times {10}^2$	$\mathbf{K}\times 10$	$\mathit{\mu }$	$\mathbf{\Sigma }$	$\mathit{\tau }$	K
5	10^3	$\underset{(2.045)}{6.646}$	$\underset{(2.673)}{12.553}$	$\underset{(19.932)}{55.538}$	$\underset{(3.540)}{12.660}$	3.249	4.696	2.786	3.576
	10^6	$\underset{(0.067)}{0.208}$	$\underset{(0.081)}{0.386}$	$\underset{(0.604)}{1.829}$	$\underset{(0.111)}{0.395}$	3.107	4.779	3.031	3.550
10	10^3	$\underset{(2.104)}{9.893}$	$\underset{(2.377)}{18.278}$	$\underset{(40.337)}{147.627}$	$\underset{(9.648)}{35.845}$	4.703	7.688	3.660	3.715
	10^6	$\underset{(0.072)}{0.313}$	$\underset{(0.079)}{0.588}$	$\underset{(1.353)}{4.860}$	$\underset{(0.224)}{1.080}$	4.370	7.480	3.591	4.826

We assume a multivariate standard normal distribution. The dimension of the system is fixed at n = 5 or 10. The second column states the sample size m in each simulation, so $m={10}^3$ or 10⁶. For each sample of size m we compute the mean, covariance, Kollo skewness, and Kollo kurtosis and record the RMSE derived using (2)(2) $$\sqrt{\sum _i^n{({\mu }_i-{\widehat{\mu }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\sigma }_{\mathit{\text{ij}}}-{\widehat{\sigma }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},\sqrt{\sum _i^n{({\tau }_i-{\widehat{\tau }}_i)}^2}\hspace{0.17em},\sqrt{\sum _{i,j}^n{({\kappa }_{\mathit{\text{ij}}}-{\widehat{\kappa }}_{\mathit{\text{ij}}})}^2}\hspace{0.17em},$$(2) . We repeat this many times and record the average RMSE and its standard error (in parentheses, below) in columns 3–6. The right-hand panel shows the ratio of the average RMSE to its standard error.

3. Conditions for matching Kollo skewness and kurtosis

Consider the simple case of generating a sample ${\mathbf{X}}_{\mathit{\text{mn}}}$ of size m from an n-dimensional distribution that has the exact $n\times 1$ mean vector ${\mathit{\mu }}_n$ and covariance matrix ${\mathbf{\Sigma }}_n.$ Let ${\mathbf{A}}_n$ be a fixed square matrix such that ${\mathbf{A}}_n{\mathbf{A}}_n^{\prime }={\mathbf{\Sigma }}_n.$ Following Ledermann et al. (2011) we can select any “scaled L-matrix” ${\mathbf{M}}_{\mathit{\text{mn}}}$ that is defined by the properties: (3) $$\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{\mathbf{1}}_m^{\prime }{\mathit{M}}_{\mathit{\text{mn}}}={\mathbf{0}}_n^{\prime }\hspace{0.17em},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{\mathit{M}}_{\mathit{\text{mn}}}^{\prime }{\mathit{M}}_{\mathit{\text{mn}}}=m{\mathit{I}}_n\hspace{0.17em},$$(3) where ${\mathbf{1}}_m$ denotes an $m\times 1$ vector of 1s, ${\mathbf{0}}_m$ is an $m\times 1$ vector of 0s and ${\mathbf{I}}_n$ is the n × n identity matrix. The matrix ${\mathbf{M}}_{\mathit{\text{mn}}}$ is just an L matrix scaled by $m^{1/2}.$ See Ledermann and Alexander (2012) for more details about different types of L matrices that can be used to target the Mardia measures exactly, either skewness or kurtosis, but not both. The original concept of ROM simulation selects any random n × n rotation matrix ${\mathbf{R}}_n,$ and possibly also a random m × m permutation matrix ${\mathbf{Q}}_m,$ so that by (3)(3) $$\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{\mathbf{1}}_m^{\prime }{\mathit{M}}_{\mathit{\text{mn}}}={\mathbf{0}}_n^{\prime }\hspace{0.17em},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{\mathit{M}}_{\mathit{\text{mn}}}^{\prime }{\mathit{M}}_{\mathit{\text{mn}}}=m{\mathit{I}}_n\hspace{0.17em},$$(3) every sample of the form: (4) $${\mathit{X}}_{\mathit{\text{mn}}}={\mathbf{1}}_m{\mathit{\mu }}_n^{\prime }+{\mathit{Q}}_m{\mathit{M}}_{\mathit{\text{mn}}}{\mathit{R}}_n{\mathit{A}}_n\hspace{0.17em},$$(4) has mean ${\mathit{\mu }}_n$ and covariance matrix ${\mathbf{\Sigma }}_n.$ The problem with such a construction is that these ROM simulations can have quite different Kollo skewness and kurtosis because they depend on the choice of ${\mathbf{R}}_n.$ For this reason, we set ${\mathbf{R}}_n={\mathbf{I}}_n$ in Equation (4)(4) $${\mathit{X}}_{\mathit{\text{mn}}}={\mathbf{1}}_m{\mathit{\mu }}_n^{\prime }+{\mathit{Q}}_m{\mathit{M}}_{\mathit{\text{mn}}}{\mathit{R}}_n{\mathit{A}}_n\hspace{0.17em},$$(4) and introduce another source of randomness.

First, we note that ${\mathbf{M}}_{\mathit{\text{mn}}}$ must satisfy more than just (3)(3) $$\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{\mathbf{1}}_m^{\prime }{\mathit{M}}_{\mathit{\text{mn}}}={\mathbf{0}}_n^{\prime }\hspace{0.17em},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{\mathit{M}}_{\mathit{\text{mn}}}^{\prime }{\mathit{M}}_{\mathit{\text{mn}}}=m{\mathit{I}}_n\hspace{0.17em},$$(3) for the simulation ${\mathbf{X}}_{\mathit{\text{mn}}}$ given by: (5) $${\mathit{X}}_{\mathit{\text{mn}}}={\mathbf{1}}_m{\mathit{\mu }}_n^{\prime }+{\mathit{Q}}_m{\mathit{M}}_{\mathit{\text{mn}}}{\mathit{A}}_n\hspace{0.17em},$$(5) to have exact Kollo skewness ${\mathit{\tau }}_n$ and Kollo kurtosis ${\mathbf{K}}_n$ as well as the exact mean vector ${\mathit{\mu }}_n$ and covariance matrix ${\mathbf{\Sigma }}_n.$ It is easiest to write these necessary and sufficient conditions in summation notation, letting ${\mathbf{m}}_i^r$ denote the $1\times n$ vector that is the ith row of ${\mathbf{M}}_{\mathit{\text{mn}}},$ as follows: (6a) $$\sum _{i=1}^m{\mathit{m}}_i^r={\mathbf{0}}_n^{\prime }\hspace{0.17em},\mathrm{\hspace{1em}}\sum _{i=1}^m{{\mathit{m}}_i^r}^{\mathbf{\prime }}{\mathit{m}}_i^r=m{\mathit{I}}_n,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\left(\text{For matching}{\mathit{\mu }}_n\text{and}{\mathbf{\Sigma }}_n\right)$$(6a) (6b) $$\sum _{i=1}^m{\left({\mathbf{1}}_n^{\prime }{{\mathit{m}}_i^r}^{\prime }\right)}^2{{\mathit{m}}_i^r}^{\prime }=m{\mathit{\tau }}_n,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\left(\text{For matching}{\mathit{\tau }}_n\right)$$(6b) (6c) $$\sum _{i=1}^m{\left({\mathbf{1}}_n^{\prime }{{\mathit{m}}_i^r}^{\prime }\right)}^2{{\mathit{m}}_i^r}^{\prime } {\mathit{m}}_i^r=m{\mathit{K}}_n.\left(\text{For matching}{\mathit{K}}_n\right)$$(6c)

Unfortunately (6b)(6b) $$\sum _{i=1}^m{\left({\mathbf{1}}_n^{\prime }{{\mathit{m}}_i^r}^{\prime }\right)}^2{{\mathit{m}}_i^r}^{\prime }=m{\mathit{\tau }}_n,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\left(\text{For matching}{\mathit{\tau }}_n\right)$$(6b) and (6c)(6c) $$\sum _{i=1}^m{\left({\mathbf{1}}_n^{\prime }{{\mathit{m}}_i^r}^{\prime }\right)}^2{{\mathit{m}}_i^r}^{\prime } {\mathit{m}}_i^r=m{\mathit{K}}_n.\left(\text{For matching}{\mathit{K}}_n\right)$$(6c) are difficult to solve directly. Instead, we follow Hanke et al. (2017) and introduce an n × n rotation matrix, ${\mathbf{\Omega }}_n,$ which satisfies ${\mathbf{1}}^{\mathbf{\prime }}{\mathbf{\Omega }}_n=\sqrt{n}[1,0,\dots ,0]$ in addition to the usual conditions ${\mathbf{\Omega }}_n^{\prime }{\mathbf{\Omega }}_n={\mathbf{I}}_n$ and $\mathrm{\text{det}}({\mathbf{\Omega }}_n)=1.$⁵ One example of such a matrix is: $${\mathbf{\Omega }}_n=\left(\begin{array}{ccccccc}\frac{1}{\sqrt{n}}& \frac{{(-1)}^{n-1}}{\sqrt{n/(n-1)}}& 0& 0& \cdots & 0& 0\\ \frac{1}{\sqrt{n}}& \frac{{(-1)}^n}{\sqrt{n(n-1)}}& \frac{{(-1)}^{n-1}}{\sqrt{(n-1)/(n-2)}}& 0& \cdots & 0& 0\\ \frac{1}{\sqrt{n}}& \frac{{(-1)}^n}{\sqrt{n(n-1)}}& \frac{{(-1)}^n}{\sqrt{(n-1)(n-2)}}& \frac{{(-1)}^{n-1}}{\sqrt{(n-2)/(n-3)}}& \cdots & 0& 0\\ \frac{1}{\sqrt{n}}& \frac{{(-1)}^n}{\sqrt{n(n-1)}}& \frac{{(-1)}^n}{\sqrt{(n-1)(n-2)}}& \frac{{(-1)}^n}{\sqrt{(n-2)(n-3)}}& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ \frac{1}{\sqrt{n}}& \frac{{(-1)}^n}{\sqrt{n(n-1)}}& \frac{{(-1)}^n}{\sqrt{(n-1)(n-2)}}& \frac{{(-1)}^n}{\sqrt{(n-2)(n-3)}}& \cdots & \frac{{(-1)}^{n-1}}{\sqrt{3/2}}& 0\\ \frac{1}{\sqrt{n}}& \frac{{(-1)}^n}{\sqrt{n(n-1)}}& \frac{{(-1)}^n}{\sqrt{(n-1)(n-2)}}& \frac{{(-1)}^n}{\sqrt{(n-2)(n-3)}}& \cdots & \frac{{(-1)}^n}{\sqrt{3\times 2}}& \frac{{(-1)}^{n-1}}{\sqrt{2}}\\ \frac{1}{\sqrt{n}}& \frac{{(-1)}^n}{\sqrt{n(n-1)}}& \frac{{(-1)}^n}{\sqrt{(n-1)(n-2)}}& \frac{{(-1)}^n}{\sqrt{(n-2)(n-3)}}& \cdots & \frac{{(-1)}^n}{\sqrt{3\times 2}}& \frac{{(-1)}^n}{\sqrt{2}}\end{array}\right).$$

The rotation matrix is used to convert ${\mathbf{M}}_{\mathit{\text{mn}}}$ to a new matrix ${\mathbf{S}}_{\mathit{\text{mn}}},$ i.e., ${\mathbf{S}}_{\mathit{\text{mn}}}={\mathbf{M}}_{\mathit{\text{mn}}}{\mathbf{\Omega }}_n,$ and so ${\mathbf{M}}_{\mathit{\text{mn}}}={\mathbf{S}}_{\mathit{\text{mn}}}{\mathbf{\Omega }}_n^{\prime }.$ The sample construction (5)(5) $${\mathit{X}}_{\mathit{\text{mn}}}={\mathbf{1}}_m{\mathit{\mu }}_n^{\prime }+{\mathit{Q}}_m{\mathit{M}}_{\mathit{\text{mn}}}{\mathit{A}}_n\hspace{0.17em},$$(5) may now be rewritten: (7) $${\mathit{X}}_{\mathit{\text{mn}}}={\mathbf{1}}_m{\mathit{\mu }}_n^{\prime }+{\mathit{Q}}_m{\mathit{S}}_{\mathit{\text{mn}}}{\mathbf{\Omega }}_n^{\prime }{\mathit{A}}_n,$$(7) and the necessary and sufficient conditions (6) become: (8a) $${\mathbf{1}}_m^{\prime }{\mathit{S}}_{\mathit{\text{mn}}}={\mathbf{0}}_n,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{\mathit{S}}_{\mathit{\text{mn}}}^{\prime }{\mathit{S}}_{\mathit{\text{mn}}}=m{\mathit{I}}_n,$$(8a) (8b) $$\sum _{i=1}^m{s}_{i1}{}^2 {{\mathit{s}}_i^r}^{\prime }=\frac{m}{n} {\mathbf{\Omega }}_n^{\prime }{\mathit{\tau }}_n:=m{\tilde{\mathit{\tau }}}_n,$$(8b) (8c) $$\sum _{i=1}^m{s}_{i1}{}^2 {{\mathit{s}}_i^r}^{\prime } {\mathit{s}}_i^r=\frac{m}{n} {\mathbf{\Omega }}_n^{\prime }{\mathit{K}}_n{\mathbf{\Omega }}_n:=m{\tilde{\mathit{K}}}_n.$$(8c)

The above formulation explains why we require ${\mathbf{1}}^{\prime }{\mathbf{\Omega }}_n=\sqrt{n}[1,0,\dots ,0].$ We can see that the most complex terms ${\left({\mathbf{1}}_n^{\prime }{{\mathbf{m}}_i^r}^{\prime }\right)}^2$ in (6b)(6b) $$\sum _{i=1}^m{\left({\mathbf{1}}_n^{\prime }{{\mathit{m}}_i^r}^{\prime }\right)}^2{{\mathit{m}}_i^r}^{\prime }=m{\mathit{\tau }}_n,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\left(\text{For matching}{\mathit{\tau }}_n\right)$$(6b) and (6c)(6c) $$\sum _{i=1}^m{\left({\mathbf{1}}_n^{\prime }{{\mathit{m}}_i^r}^{\prime }\right)}^2{{\mathit{m}}_i^r}^{\prime } {\mathit{m}}_i^r=m{\mathit{K}}_n.\left(\text{For matching}{\mathit{K}}_n\right)$$(6c) are reduced to terms of the form $n{s}_{i1}^2$ in (8b)(8b) $$\sum _{i=1}^m{s}_{i1}{}^2 {{\mathit{s}}_i^r}^{\prime }=\frac{m}{n} {\mathbf{\Omega }}_n^{\prime }{\mathit{\tau }}_n:=m{\tilde{\mathit{\tau }}}_n,$$(8b) and (8c)(8c) $$\sum _{i=1}^m{s}_{i1}{}^2 {{\mathit{s}}_i^r}^{\prime } {\mathit{s}}_i^r=\frac{m}{n} {\mathbf{\Omega }}_n^{\prime }{\mathit{K}}_n{\mathbf{\Omega }}_n:=m{\tilde{\mathit{K}}}_n.$$(8c) , with $s_{i1}$ denoting the first element in the ith row of ${\mathbf{S}}_{\mathit{\text{mn}}}.$ This way, the focus moves from simulating (5)(5) $${\mathit{X}}_{\mathit{\text{mn}}}={\mathbf{1}}_m{\mathit{\mu }}_n^{\prime }+{\mathit{Q}}_m{\mathit{M}}_{\mathit{\text{mn}}}{\mathit{A}}_n\hspace{0.17em},$$(5) using different ${\mathbf{M}}_n$ to simulating (7)(7) $${\mathit{X}}_{\mathit{\text{mn}}}={\mathbf{1}}_m{\mathit{\mu }}_n^{\prime }+{\mathit{Q}}_m{\mathit{S}}_{\mathit{\text{mn}}}{\mathbf{\Omega }}_n^{\prime }{\mathit{A}}_n,$$(7) using different ${\mathbf{S}}_{\mathit{\text{mn}}}$ which satisfies (8).

4. Theoretical results

For convenience, we now rewrite the conditions (8) as functions of each column of ${\mathbf{S}}_{\mathit{\text{mn}}},$ respectively, using the Hadamard product $°.$ Let ${\mathbf{s}}_j^c$ denote the $m\times 1$ vector that is the jth column of ${\mathbf{S}}_{\mathit{\text{mn}}}.$ The following formulation of the necessary and sufficient conditions for the simulated sample to match exact mean, covariance and Kollo skewness and kurtosis will be useful in this section: (9a) $${\mathbf{1}}_m^{\prime }{\mathit{s}}_k^c=0\hspace{0.17em},$$(9a) (9b) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_j^c°{\mathit{s}}_k^c)=0\hspace{0.17em},\forall j<k\hspace{0.17em},$$(9b) (9c) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_k^c°{\mathit{s}}_k^c)=m\hspace{0.17em},$$(9c) (9d) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_k^c)=m{\tilde{\tau }}_k\hspace{0.17em},$$(9d) (9e) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c°{\mathit{s}}_k^c)=m{\tilde{\kappa }}_{\mathit{\text{ik}}},\forall i\le k,$$(9e) where ${\tilde{\tau }}_k$ is the kth element of the scaled and rotated Kollo skewness vector ${\tilde{\mathbf{\tau }}}_{\mathbf{n}}$ and ${\tilde{\kappa }}_{\mathit{\text{ik}}}$ is the element in the ith row and kth column of the scaled and rotated Kollo kurtosis matrix ${\tilde{\mathbf{K}}}_n.$ Matching the exact mean, covariance and Kollo skewness relies on solving conditions (9a)(9a) $${\mathbf{1}}_m^{\prime }{\mathit{s}}_k^c=0\hspace{0.17em},$$(9a) –(9d)(9d) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_k^c)=m{\tilde{\tau }}_k\hspace{0.17em},$$(9d) for ${\mathbf{S}}_{\mathit{\text{mn}}}.$ The additional constraints (9e)(9e) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c°{\mathit{s}}_k^c)=m{\tilde{\kappa }}_{\mathit{\text{ik}}},\forall i\le k,$$(9e) set out the k restrictions on the kth column of ${\mathbf{S}}_{\mathit{\text{mn}}}$ that must hold in order that each simulation ${\mathbf{X}}_{\mathit{\text{mn}}}$ matches the target Kollo kurtosis.

The system (9) is not easy to solve. When k = 1 we have to solve a system of linear equations for the mean constraint (9a)(9a) $${\mathbf{1}}_m^{\prime }{\mathit{s}}_k^c=0\hspace{0.17em},$$(9a) , a system of quadratic equations for the variance constraint (9c)(9c) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_k^c°{\mathit{s}}_k^c)=m\hspace{0.17em},$$(9c) , a system of cubic equations for the Kollo skewness constraint (9d)(9d) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_k^c)=m{\tilde{\tau }}_k\hspace{0.17em},$$(9d) , and a system of quartic equations for the Kollo kurtosis constraint (9e)(9e) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c°{\mathit{s}}_k^c)=m{\tilde{\kappa }}_{\mathit{\text{ik}}},\forall i\le k,$$(9e) . When k > 1 we have four different systems of linear equations, for the mean constraint (9a)(9a) $${\mathbf{1}}_m^{\prime }{\mathit{s}}_k^c=0\hspace{0.17em},$$(9a) , the covariance constraint (9b)(9b) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_j^c°{\mathit{s}}_k^c)=0\hspace{0.17em},\forall j<k\hspace{0.17em},$$(9b) , the Kollo skewness constraint (9d)(9d) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_k^c)=m{\tilde{\tau }}_k\hspace{0.17em},$$(9d) , and Kollo kurtosis constraints (9e)(9e) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c°{\mathit{s}}_k^c)=m{\tilde{\kappa }}_{\mathit{\text{ik}}},\forall i\le k,$$(9e) when i < k, and two different systems of quadratic equations for the variance constraint (9c)(9c) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_k^c°{\mathit{s}}_k^c)=m\hspace{0.17em},$$(9c) and the Kollo kurtosis constraint (9e)(9e) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c°{\mathit{s}}_k^c)=m{\tilde{\kappa }}_{\mathit{\text{ik}}},\forall i\le k,$$(9e) when i = k. In total, there are $n(n+5)/2$ equations to solve for matching the first three moments, or $n(n+3)$ equations for matching the first four moments. For example, in a 10-dimensional system, we need to solve 75 equations for matching the first three moments, or 130 equations to match the first four moments.

In this section we show how to solve (9) for ${\mathbf{S}}_{\mathit{\text{mn}}}$ column by column, by (a) finding a particular type of orthonormal basis for the sample space of the n random variables; and then (b) solving for the coefficients which generate a subspace ${\mathbf{S}}_{\mathit{\text{mn}}}$ of its span which contains samples having exactly the target Kollo moments. That is, for each column k and for some column dimension l_k, with $1\le l_k\le m,$ we first find a specific $m\times l_k$ rectangular orthonormal matrix ${\mathbf{\Theta }}_{m,l_k}^k,$ i.e., a matrix having the properties $${\mathbf{1}}_m^{\prime }{\mathbf{\Theta }}_{m,l_k}^k={\mathbf{0}}_{l_k}\hspace{1em}\text{and}\hspace{1em}{\left({\mathbf{\Theta }}_{m,l_k}^k\right)}^{\prime }{\mathbf{\Theta }}_{m,l_k}^k={\mathit{I}}_{l_k}.$$

Then we set (10) $${\mathit{s}}_k^c={\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k,$$(10) where ${\mathbf{c}}_k$ is an $l_k\times 1$ column vector with at least one non-zero element. The decomposition (10)(10) $${\mathit{s}}_k^c={\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k,$$(10) allows one to select the orthonormal matrix ${\mathbf{\Theta }}_{m,l_k}^k$ in such a way that ${\mathbf{s}}_k^c$ always satisfies Equations (9a)(9a) $${\mathbf{1}}_m^{\prime }{\mathit{s}}_k^c=0\hspace{0.17em},$$(9a) and (9b)(9b) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_j^c°{\mathit{s}}_k^c)=0\hspace{0.17em},\forall j<k\hspace{0.17em},$$(9b) . That is, instead of solving (9) for ${\mathbf{s}}_k^c$ directly, we may solve for ${\mathbf{\Theta }}_{m,l_k}^k$ and ${\mathbf{c}}_k$ iteratively, for $k=1,\dots ,n,$ using the following result:

Lemma 1.

Suppose the first k − 1 columns in the matrix ${\mathbf{S}}_{\mathit{\text{mn}}}$ satisfy conditions (9a)(9a) $${\mathbf{1}}_m^{\prime }{\mathit{s}}_k^c=0\hspace{0.17em},$$(9a) and (9b)(9b) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_j^c°{\mathit{s}}_k^c)=0\hspace{0.17em},\forall j<k\hspace{0.17em},$$(9b) and the kth column ${\mathbf{s}}_k^c$ is constructed using (10)(10) $${\mathit{s}}_k^c={\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k,$$(10) . Then the conditions (9a)(9a) $${\mathbf{1}}_m^{\prime }{\mathit{s}}_k^c=0\hspace{0.17em},$$(9a) and (9b)(9b) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_j^c°{\mathit{s}}_k^c)=0\hspace{0.17em},\forall j<k\hspace{0.17em},$$(9b) will be satisfied if the column vectors in ${\mathbf{\Theta }}_{m,l_k}^k$ are also orthogonal to the vectors ${\mathbf{1}}_m$ and ${\mathbf{s}}_1^c,\dots ,{\mathbf{s}}_{k-1}^c.$

Generating ${\mathbf{S}}_{\mathit{\text{mn}}}$ column by column this way considerably reduces the number of constraints that must be satisfied. Later we shall show that, in addition to the mean and covariance, we can match the Kollo skewness and kurtosis by placing only k + 2 constraints on the kth column of ${\mathbf{S}}_{\mathit{\text{mn}}}.$ And if we only wish to target the first three moments, there are just two constraints on each column of ${\mathbf{S}}_{\mathit{\text{mn}}}.$ In total, instead of $n(n+5)/2$ equations to solve for matching the first three moments, or $n(n+3)$ equations for the first four moments, now we only need to solve 2n and $n(n+5)/2$ equations, respectively. For instance, in a 10-dimensional system, we need to solve 20 equations to match the mean, covariance, and Kollo skewness or 75 equations to additionally match the Kollo kurtosis. Of course, the dimension reduction increases with n; e.g., when n = 20 the first three moments constraints reduce from 250 to 40, and the first four moments constraints reduce from 460 to 250.

Our method of solution requires that we fix the first k − 1 columns of ${\mathbf{S}}_{\mathit{\text{mn}}}$ before we generate the matrix ${\mathbf{\Theta }}_{m,l_k}^k.$ We then generate a rectangular orthonormal matrix with dimensions m × N whose first k columns are equal to $m^{-1/2}{\mathbf{1}}_m$ and $m^{-1/2}{\mathbf{s}}_1^c,\dots ,m^{-1/2}{\mathbf{s}}_{k-1}^c.$ The columns of the matrix ${\mathbf{\Theta }}_{m,l_k}^k$ are then randomly selected from the other N—k columns of this matrix. That is, we take a random sample (without replacement) of any column j for $j\ge k$ and the selected columns ${\mathit{\theta }}_j$ for $j=1,\dots ,l_k$ are then used to form the orthonormal matrix ${\mathbf{\Theta }}_{m,l_k}^k.$ Note that we require $N\le m,$ because there are m elements in each column and so, no more than m columns can be mutually orthogonal. This places an upper bound on l_k because $l_k\le N-k\le m-k.$

There is a close relationship between ${\mathbf{\Theta }}_{m,l_k}^k$ and the L matrices (Ledermann et al., 2011). In fact, we can immediately obtain ${\mathbf{\Theta }}_{m,l_k}^k$ by randomly extracting l_k columns from the last $N-(k-1)$ columns of an L matrix whose first k − 1 columns are $m^{-1/2}{\mathbf{s}}_1^c,\dots ,m^{-1/2}{\mathbf{s}}_{k-1}^c.$

Now we turn to the actual solution of (9) using the construction (10)(10) $${\mathit{s}}_k^c={\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k,$$(10) , starting with the simpler case of matching the first three moments only. Recall that ${\mathbf{s}}_1^c$ is different from the other columns because it involves a cubic equation, which is easy to see by rewriting the conditions (9a)—(9d) for ${\mathbf{s}}_1^c$ as: (11) $$\sum _{i=1}^m{s}_{i1}=0\hspace{0.17em},\hspace{1em}\sum _{i=1}^m{s}_{i1}{}^2=m\hspace{0.17em},\hspace{1em}\sum _{i=1}^m{s}_{i1}{}^3=m{\tilde{\tau }}_1.$$(11)

The following Theorem explains how to solve for ${\mathbf{s}}_1^c,$ the first column of ${\mathbf{S}}_{\mathit{\text{mn}}},$ from which the other columns are generated in an iterative manner:

Theorem 1.

Suppose the kth column of ${\mathbf{S}}_{\mathit{\text{mn}}}$ satisfies (10)(10) $${\mathit{s}}_k^c={\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k,$$(10) for $l_k\ge 2$, that ${\mathbf{\Theta }}_{m,l_k}^k=\left[{\mathit{\theta }}_1^k,\dots ,{\mathit{\theta }}_{l_k}^k\right]$ satisfies the conditions of Lemma 1, and that ${\mathbf{s}}_k^1,\dots ,{\mathbf{s}}_{k-1}^c$ satisfy conditions (9a)(9a) $${\mathbf{1}}_m^{\prime }{\mathit{s}}_k^c=0\hspace{0.17em},$$(9a) to (9d)(9d) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_k^c)=m{\tilde{\tau }}_k\hspace{0.17em},$$(9d) . Then: by Lemma 1, ${\mathbf{s}}_k^c$ satisfies conditions (9a)(9a) $${\mathbf{1}}_m^{\prime }{\mathit{s}}_k^c=0\hspace{0.17em},$$(9a) and (9b)(9b) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_j^c°{\mathit{s}}_k^c)=0\hspace{0.17em},\forall j<k\hspace{0.17em},$$(9b) ; it also satisfies condition (9c)(9c) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_k^c°{\mathit{s}}_k^c)=m\hspace{0.17em},$$(9c) if(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) and the condition (9d)(9d) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_k^c)=m{\tilde{\tau }}_k\hspace{0.17em},$$(9d) also holds if(13a) $$\text{for}\hspace{1em}k=1:{\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)]=m{\tilde{\tau }}_1\hspace{0.17em},$$(13a) (13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b) where$${\mathit{B}}_{m,l_k}^k=[{\mathit{b}}_1^k,\dots .{\mathit{b}}_{l_k}^k]=[{\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{\theta }}_1^k,\dots ,{\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{\theta }}_{l_k}^k].$$

The pseudocode in algorithm 1 summarizes how to generate a sample ${\mathbf{X}}_{\mathit{\text{mn}}}$ that matches any given mean vector, covariance matrix, and Kollo skewness vector.

Algorithm 1.

Matching the First Three Moments

1: function ROM₃($m,\mathit{\mu },\mathbf{\Sigma },\mathit{\tau }$)      ⊳ m: Number of observations, $\mathit{\mu },\mathbf{\Sigma },\mathit{\tau }:$ Targeted first three moments

2:      $n=\mathit{\text{length}}(\mathit{\tau })$

3:      Obtain ${\mathbf{A}}_n$ where ${\mathbf{A}}_n{\mathbf{A}}_n^{\prime }=\mathbf{\Sigma }$

4:      Generate ${\mathbf{\Omega }}_n$            ⊳ ${\mathbf{\Omega }}_n$ is an n × n rotation matrix satisfying ${\mathbf{1}}^{\prime }{\mathbf{\Omega }}_n=\sqrt{n}[1,0,\dots ,0]$

5:      $\tilde{\mathit{\tau }}={\mathbf{\Omega }}_n^{-1}{n}^{-1}\mathit{\tau }$                      ⊳ Obtain the rotated Kollo skewness

6:    while Theorem 1 conditions (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) and (13a)(13a) $$\text{for}\hspace{1em}k=1:{\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)]=m{\tilde{\tau }}_1\hspace{0.17em},$$(13a) for ${\mathbf{c}}_1$ do not hold do    ⊳ Solve for the first column

7:     Generate l₁ orthonormal vectors to form ${\mathbf{\Theta }}_{m,l_1}^1$ satisfying Lemma 1

8:     Obtain real solutions for the coefficient ${\mathbf{c}}_1$ restricted by Equations (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) and (13a)(13a) $$\text{for}\hspace{1em}k=1:{\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)]=m{\tilde{\tau }}_1\hspace{0.17em},$$(13a)

9:    Obtain ${\mathbf{s}}_1^c={\mathbf{\Theta }}_{m,l_1}^1{\mathbf{c}}_1$

10:     for $k=2,\dots ,n$ do:

11:      while Theorem 1 conditions (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) and (13b)(13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b) for ${\mathbf{c}}_k$ do not hold do

12:        Generate l_k orthonormal vectors to form ${\mathbf{\Theta }}_{m,l_k}^k$ satisfying Lemma 1

13:        Obtain real solutions for the coefficients ${\mathbf{c}}_k$ restricted by Equations (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) and (13b)(13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b)

14:      Obtain ${\mathbf{s}}_k^c={\mathbf{\Theta }}_{m,l_k}^k{\mathbf{c}}_k$

15:     Generate permutation matrix ${\mathbf{Q}}_m$ if required

16:     ${\mathbf{X}}_{\mathit{\text{mn}}}={\mathbf{1}}_m{\mathit{\mu }}^{\prime }+{\mathbf{Q}}_m{\mathbf{S}}_{\mathit{\text{mn}}}{{\Omega }^{\prime }}_n{\mathbf{A}}_n$        ⊳ Generate a sample with exact first three moments

17:     return ${\mathbf{X}}_{\mathit{\text{mn}}}$

It solves ${\mathbf{S}}_{\mathit{\text{mn}}}$ iteratively, column-by-column while building the columns of an orthonormal matrix ${\mathbf{\Theta }}_{m,l_k}^k$ that yields a real-valued solution for ${\mathbf{c}}_k.$ This means, to obtain ${\mathbf{s}}_k^c$ satisfying conditions (9a)—(9d), the k − 1 columns of ${\mathbf{S}}_{\mathit{\text{mn}}}$ must first be known and will satisfy conditions (9a)—(9d). So we first obtain an orthonormal matrix ${\mathbf{\Theta }}_{m,l_1}^1$ satisfying Lemma 1 and then obtain the coefficients ${\mathbf{c}}_1$ by solving a system of quadratic Equation (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) and a system of cubic Equation (13a)(13a) $$\text{for}\hspace{1em}k=1:{\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)]=m{\tilde{\tau }}_1\hspace{0.17em},$$(13a) simultaneously. Next, we generate ${\mathbf{s}}_k^c$ for $k=2,\dots ,n.$ Now ${\mathbf{c}}_k$ should satisfy a quadratic system (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) and a linear system (13b)(13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b) simultaneously.

We require $l_k\ge 2$ because the elements of ${\mathbf{c}}_k$ are always restricted by 2 equations: for k = 1 the systems (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) , (13a)(13a) $$\text{for}\hspace{1em}k=1:{\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)]=m{\tilde{\tau }}_1\hspace{0.17em},$$(13a) —and for k > 1 the systems (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) and (13b)(13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b) —are both exactly determined when $l_k=2;$ but when $l_k>2$ the systems are under-determined, which makes it easier to find a real-value solution for ${\mathbf{c}}_k.$

While both general and flexible, there is an element of trial and error for finding a real-valued solution in Algorithm 1. The following Corollary 1 considers a special case of Theorem 1 for l_k = 2 and $k\ge 2.$ It is useful because it derives a single criterion for the existence of a real-valued solution to the system of (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) and (13b)(13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b) and an analytical formula for each 2 × 1 vector ${\mathbf{c}}_k={\left[c_1^k,c_2^k\right]}^{\prime }$ when such a solution exists.

Corollary 1.

When l_k = 2 and $k\ge 2,$ the system of Equations (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) and (13b)(13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b) has a real solution if (14) $${\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,2}^k {{\mathit{B}}_{m,2}^k}^{\prime }{\mathbf{1}}_m\ge m{{\tilde{\tau }}_k}^2\hspace{0.17em},$$(14) and when the condition (14)(14) $${\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,2}^k {{\mathit{B}}_{m,2}^k}^{\prime }{\mathbf{1}}_m\ge m{{\tilde{\tau }}_k}^2\hspace{0.17em},$$(14) holds the solution is:(15) $$c_1^k=\frac{{\tilde{b}}_1^{k}m{\tilde{\tau }}_k\pm {\tilde{b}}_2^k\sqrt{\left[{\left({\tilde{b}}_1^k\right)}^2\mathrm{}+\mathrm{}{\left({\tilde{b}}_2^k\right)}^2\right]m\mathrm{}‐\mathrm{}{m}^2{\tilde{\tau }}_k{}^2}}{{\left({\tilde{b}}_1^k\right)}^2\mathrm{}+\mathrm{}{\left({\tilde{b}}_2^k\right)}^2}\hspace{0.17em},\hspace{1em}{c}_2^k=\frac{{\tilde{b}}_2^{k}m{\tilde{\tau }}_k\mp {\tilde{b}}_1^k\sqrt{\left[{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2\right]m\mathrm{}‐\mathrm{}{m}^2{\tilde{\tau }}_k{}^2}}{{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2}\hspace{0.17em},$$(15) where ${\tilde{b}}_1^k$ and ${\tilde{b}}_2^k$ are the first and second elements of ${\mathbf{1}}_m^{\prime }{\mathbf{B}}_{m,2}^k.$

Corollary 1 greatly improves the computational efficiency of the new method. Setting l_k = 2 for k > 1 means that the matrix has only two columns which are simple to orthogonalize. Moreover, it is very quick to repeatedly generate orthonormal matrices ${\mathbf{\Theta }}_{m,2}^k$ for k > 1 until we find those which satisfy the criterion (14)(14) $${\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,2}^k {{\mathit{B}}_{m,2}^k}^{\prime }{\mathbf{1}}_m\ge m{{\tilde{\tau }}_k}^2\hspace{0.17em},$$(14) , and then we can immediately obtain values for ${\mathbf{c}}_k$ using the analytical formula (15)(15) $$c_1^k=\frac{{\tilde{b}}_1^{k}m{\tilde{\tau }}_k\pm {\tilde{b}}_2^k\sqrt{\left[{\left({\tilde{b}}_1^k\right)}^2\mathrm{}+\mathrm{}{\left({\tilde{b}}_2^k\right)}^2\right]m\mathrm{}‐\mathrm{}{m}^2{\tilde{\tau }}_k{}^2}}{{\left({\tilde{b}}_1^k\right)}^2\mathrm{}+\mathrm{}{\left({\tilde{b}}_2^k\right)}^2}\hspace{0.17em},\hspace{1em}{c}_2^k=\frac{{\tilde{b}}_2^{k}m{\tilde{\tau }}_k\mp {\tilde{b}}_1^k\sqrt{\left[{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2\right]m\mathrm{}‐\mathrm{}{m}^2{\tilde{\tau }}_k{}^2}}{{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2}\hspace{0.17em},$$(15) . In this way, we do not need to search for the numerical solution for the simultaneous Equations (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) and (13b)(13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b) .

Next, we consider simulations matching mean, covariance, Kollo skewness, and Kollo kurtosis. When matching only the first three moments, there is no reason why we should not always set l_k = 2 and use Corollary 1. Then the moment matching algorithm is very fast, being fully analytic. However, the next Theorem, which explains how to match the first four moments, requires $l_k\ge k+2.$ Again we begin by generating the first column of ${\mathbf{S}}_{\mathit{\text{mn}}}$ because it requires higher-order conditions than the subsequent columns. To target the first four moments ${\mathbf{s}}_1^c$ should satisfy not only conditions (11)(11) $$\sum _{i=1}^m{s}_{i1}=0\hspace{0.17em},\hspace{1em}\sum _{i=1}^m{s}_{i1}{}^2=m\hspace{0.17em},\hspace{1em}\sum _{i=1}^m{s}_{i1}{}^3=m{\tilde{\tau }}_1.$$(11) but also a quartic condition for matching Kollo kurtosis: (16) $$\sum _{i=1}^m{s}_{i1}{}^4=m{\tilde{\kappa }}_{11}\hspace{0.17em},$$(16) where ${\tilde{\kappa }}_{11}$ is the first element in the first row of ${\tilde{\mathbf{K}}}_n.$ Theorem 2 provides the necessary and sufficient conditions on ${\mathbf{c}}_k$ for ${\mathbf{S}}_{\mathit{\text{mn}}}$ to satisfy conditions (9).

Theorem 2.

Suppose the kth column of ${\mathbf{S}}_{\mathit{\text{mn}}}$ satisfies (10)(10) $${\mathit{s}}_k^c={\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k,$$(10) for $l_k\ge k+2$, that ${\mathbf{\Theta }}_{m,l_k}^k=\left[{\mathit{\theta }}_1^k,\dots ,{\mathit{\theta }}_{l_k}^k\right]$ satisfies the conditions of Lemma 1, and that ${\mathbf{s}}_k^1,\dots ,{\mathbf{s}}_{k-1}^c$ satisfy conditions (9). Then ${\mathbf{s}}_k^c$ satisfies:

1. the conditions (9a)(9a) $${\mathbf{1}}_m^{\prime }{\mathit{s}}_k^c=0\hspace{0.17em},$$(9a) and (9b)(9b) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_j^c°{\mathit{s}}_k^c)=0\hspace{0.17em},\forall j<k\hspace{0.17em},$$(9b) under Lemma 1;

2. the condition (9c)(9c) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_k^c°{\mathit{s}}_k^c)=m\hspace{0.17em},$$(9c) if condition (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) holds;

3. the condition (9d)(9d) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_k^c)=m{\tilde{\tau }}_k\hspace{0.17em},$$(9d) if condition (13) holds;

4. the condition (9e)(9e) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c°{\mathit{s}}_k^c)=m{\tilde{\kappa }}_{\mathit{\text{ik}}},\forall i\le k,$$(9e) if

for k = 1: (17a) $${\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)]=m{\tilde{\kappa }}_{11}\hspace{0.17em},$$(17a) for $k=2,\dots ,n$:(17b) $${\mathit{c}}_k^{\prime }{{\mathit{D}}_{m,l_k}^k}^{\prime }{\mathit{D}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\kappa }}_{\mathit{\text{kk}}}\hspace{0.17em},$$(17b) (17c) $$\text{and}\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{{\mathit{H}}_{m,k-1}^k}^{\prime }{\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\mathit{\kappa }}}_{k-1},$$(17c) where the $(k-1)\times 1$ column vector ${\tilde{\mathit{\kappa }}}_{k-1}$ contains the first k − 1 elements of the kth column of ${\tilde{\mathbf{K}}}_n,$ ${\mathbf{D}}_{m,l_k}^k=[{\mathbf{s}}_1^c°{\mathit{\theta }}_1^k\hspace{0.17em},\dots ,{\mathbf{s}}_1^c°{\mathit{\theta }}_{l_k}^k]$ and ${\mathbf{H}}_{m,k-1}^k=[{\mathbf{s}}_1^c°{\mathbf{s}}_1^c°{\mathbf{s}}_1^c\hspace{0.17em},\dots ,{\mathbf{s}}_1^c°{\mathbf{s}}_1^c°{\mathbf{s}}_1^{k-1}].$

Algorithm 2 is the pseudocode for implementing the above results. It yields a realisation ${\mathbf{X}}_{\mathit{\text{mn}}}$ of m observations on n random variables which matches any given target mean, covariance, Kollo skewness, and Kollo kurtosis. To obtain ${\mathbf{s}}_k^c$ satisfying all the conditions (9), the k − 1 columns of ${\mathbf{S}}_{\mathit{\text{mn}}}$ should be known and have satisfied condition (9). For this, we first obtain an orthonormal matrix ${\mathbf{\Theta }}_{m,l_1}^1$ satisfying Lemma 1 and then obtain the coefficients ${\mathbf{c}}_1$ by solving a quadratic (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) , a cubic (13a)(13a) $$\text{for}\hspace{1em}k=1:{\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)]=m{\tilde{\tau }}_1\hspace{0.17em},$$(13a) , and a quartic (17a)(17a) $${\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)]=m{\tilde{\kappa }}_{11}\hspace{0.17em},$$(17a) simultaneously. Next, we generate ${\mathbf{s}}_k^c$ in an iterative manner from k = 2 to k = n. Now the vector ${\mathbf{c}}_k$ should satisfy k + 2 equations including two quadratic Equations (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) and (17b)(17b) $${\mathit{c}}_k^{\prime }{{\mathit{D}}_{m,l_k}^k}^{\prime }{\mathit{D}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\kappa }}_{\mathit{\text{kk}}}\hspace{0.17em},$$(17b) and k linear Equations (13b)(13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b) and (17c)(17c) $$\text{and}\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{{\mathit{H}}_{m,k-1}^k}^{\prime }{\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\mathit{\kappa }}}_{k-1},$$(17c) , simultaneously.

Algorithm 2.

Matching the First Four Moments

1: function ROM₄($m,\mathit{\mu },\mathbf{\Sigma },\mathit{\tau },\mathbf{K}$)    ⊳ m: Number of observations, $\mathit{\mu },\mathbf{\Sigma },\mathit{\tau },\mathbf{K}:$ Targeted first four moments

2:   $n=\mathit{\text{length}}(\mathit{\mu })$

3:   Obtain ${\mathbf{A}}_n$ where ${\mathbf{A}}_n{\mathbf{A}}_n^{\prime }=\mathbf{\Sigma }$

4:   Generate ${\mathbf{\Omega }}_n$           ⊳ ${\mathbf{\Omega }}_n$ is an n × n rotation matrix satisfying $\mathbf{1}\prime {\mathbf{\Omega }}_n=\sqrt{n}\left[1,0,\dots ,0\right]$

5:   $\tilde{\mathit{\tau }}={\mathbf{\Omega }}_n^{-1}{n}^{-1}\mathit{\tau }$                      ⊳ Obtain the rotated Kollo skewness

6:   $\tilde{\mathbf{K}}={\mathbf{\Omega }}_n^{-1}{n}^{-1}\mathbf{K}{\mathbf{\Omega }}_n$                     ⊳ Obtain the rotated Kollo kurtosis

7:   while Theorem 2 does not hold do                 ⊳ Solve for the first column

8:     generate l₁ orthonormal vectors and form ${\mathbf{\Theta }}_{m,l_1}^1$ satisfying Lemma 1

9:     obtain real solutions for the coefficient ${\mathbf{c}}_1$ restricted by Equations (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) , (13a)(13a) $$\text{for}\hspace{1em}k=1:{\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)]=m{\tilde{\tau }}_1\hspace{0.17em},$$(13a) , and (17a)(17a) $${\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)]=m{\tilde{\kappa }}_{11}\hspace{0.17em},$$(17a)

10:   Obtain ${\mathbf{s}}_1^c={\mathbf{\Theta }}_{m,l_1}^1{\mathbf{c}}_1$

11:   for $k=2,\dots ,n$ do ⊳ Solve for the other columns

12:    while Theorem 2 does not hold do

13:      Generate l_k orthonormal vectors to form ${\mathbf{\Theta }}_{m,l_k}^k$ satisfying Lemma 1

14:      Obtain real solutions for the coefficients ${\mathbf{c}}_k$ restricted by Equations (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) , (13b)(13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b) , (17b)(17b) $${\mathit{c}}_k^{\prime }{{\mathit{D}}_{m,l_k}^k}^{\prime }{\mathit{D}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\kappa }}_{\mathit{\text{kk}}}\hspace{0.17em},$$(17b) , and (17c)(17c) $$\text{and}\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{{\mathit{H}}_{m,k-1}^k}^{\prime }{\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\mathit{\kappa }}}_{k-1},$$(17c)

15:    Obtain ${\mathbf{s}}_k^c={\mathbf{\Theta }}_{m,l_k}^k{\mathbf{c}}_k$

16:   Generate permutation matrix ${\mathbf{Q}}_m$ if required

17:   ${\mathbf{X}}_{\mathit{\text{mn}}}={\mathbf{1}}_m\mathit{\mu }\prime +{\mathbf{Q}}_m{\mathbf{S}}_{\mathit{\text{mn}}}\Omega {\prime }_n{\mathbf{A}}_n$ ⊳ Generate sample with exact first four moments

18:   return ${\mathbf{X}}_{\mathit{\text{mn}}}$

We require $l_k\ge k+2$ because the elements of ${\mathbf{c}}_k$ are always restricted by k + 2 equations: for k = 1 the systems (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) , (13a)(13a) $$\text{for}\hspace{1em}k=1:{\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)]=m{\tilde{\tau }}_1\hspace{0.17em},$$(13a) , (17a)(17a) $${\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)]=m{\tilde{\kappa }}_{11}\hspace{0.17em},$$(17a) —and for k > 1 the systems (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) , (13b)(13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b) , (17b)(17b) $${\mathit{c}}_k^{\prime }{{\mathit{D}}_{m,l_k}^k}^{\prime }{\mathit{D}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\kappa }}_{\mathit{\text{kk}}}\hspace{0.17em},$$(17b) , (17c)(17c) $$\text{and}\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{{\mathit{H}}_{m,k-1}^k}^{\prime }{\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\mathit{\kappa }}}_{k-1},$$(17c) —are both exactly determined when $l_k=k+2;$ but when $l_k>k+2$ the systems are under-determined, which makes it easier to find a real-value solution for ${\mathbf{c}}_k.$

Finally, we note the following property for Kollo kurtosis:⁶ Let ${\mathbf{M}}^{(k)}$ be the matrix of dimension $m_k\times n$ with $m_k\ge n+2$ satisfying conditions (3)(3) $$\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{\mathbf{1}}_m^{\prime }{\mathit{M}}_{\mathit{\text{mn}}}={\mathbf{0}}_n^{\prime }\hspace{0.17em},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{\mathit{M}}_{\mathit{\text{mn}}}^{\prime }{\mathit{M}}_{\mathit{\text{mn}}}=m{\mathit{I}}_n\hspace{0.17em},$$(3) . Denote the ith row of ${\mathbf{M}}^{(k)}$ by ${\mathbf{m}}_i^{(k)r},$ for $i=1,\dots ,N$ where N is the number of concatenations. Then the Kollo kurtosis of the concatenated sample is the weighted sum of their Kollo kurtoses: (18) $$\begin{array}{cc}\mathit{K}\left(\left[\begin{array}{c}{\mathit{M}}^{\left(1\right)}\\ ⋮\\ {\mathit{M}}^{\left(N\right)}\end{array}\right]\right)=& \frac{\sum _{i=1}^{m_1}{\left({\mathit{m}}_i^{\left(1\right)r}{\mathbf{1}}_n\right)}^2{{\mathit{m}}_i^{\left(1\right)r}}^{\prime }{\mathit{m}}_i^{\left(1\right)r}\mathrm{}+\mathrm{}\dots \mathrm{}+\mathrm{}\sum _{i=1}^{m_N}{\left({\mathit{m}}_i^{\left(N\right)r}{\mathbf{1}}_n\right)}^2{{\mathit{m}}_i^{\left(N\right)r}}^{\prime }{\mathit{m}}_i^{\left(N\right)r}}{m_1\mathrm{}+\mathrm{}\dots \mathrm{}+\mathrm{}{m}_N}\\ =& \frac{m_1\mathit{K}\left({\mathit{M}}^{\left(1\right)}\right)+\dots +m_N\mathit{K}\left({\mathit{M}}^{\left(N\right)}\right)}{m_1\mathrm{}+\mathrm{}\dots \mathrm{}+\mathrm{}{m}_N}.\end{array}$$(18)

Hence, Kollo kurtosis is invariant under sample concatenation if each sub-matrix ${\mathbf{M}}^{(k)}$ has identical Kollo kurtosis.

5. Discussions

This paper derives necessary and sufficient analytic conditions for exactly matching the mean, covariance, Kollo skewness, and Kollo kurtosis of multivariate random variables. The necessary and sufficient conditions (9) are a highly under-determined system. Even without the condition for matching Kollo kurtosis (9e)(9e) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c°{\mathit{s}}_k^c)=m{\tilde{\kappa }}_{\mathit{\text{ik}}},\forall i\le k,$$(9e) , it is already a complex problem to find columns of ${\mathbf{S}}_{\mathit{\text{mn}}}$ which satisfy (9a)(9a) $${\mathbf{1}}_m^{\prime }{\mathit{s}}_k^c=0\hspace{0.17em},$$(9a) —(9d) because of the skewness condition (9d)(9d) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_k^c)=m{\tilde{\tau }}_k\hspace{0.17em},$$(9d) . This explains why Kollo skewness matching methods, such as the algorithm proposed by Alexander et al. (2022), are much slower than those matching only the first two moments exactly, such as Ledermann et al. (2011). Adding the conditions (9e)(9e) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c°{\mathit{s}}_k^c)=m{\tilde{\kappa }}_{\mathit{\text{ik}}},\forall i\le k,$$(9e) creates even more of a challenge: there are already k + 2 constraints on the kth column for Kollo skewness matching and an extra set of k conditions for Kollo kurtosis matching, so that in total for each column k, for $k=1,\dots ,n,$ there are $2k+2$ constraints on m observations making a total of $2n(n+1)/2+2n=n(n+3)$ conditions.

In response to this issue, we propose the novel idea of decomposing ${\mathbf{S}}_{\mathit{\text{mn}}}$ in (7)(7) $${\mathit{X}}_{\mathit{\text{mn}}}={\mathbf{1}}_m{\mathit{\mu }}_n^{\prime }+{\mathit{Q}}_m{\mathit{S}}_{\mathit{\text{mn}}}{\mathbf{\Omega }}_n^{\prime }{\mathit{A}}_n,$$(7) into the product of an orthonormal matrix and a vector of coefficients. We use the necessary and sufficient conditions to solve for the coefficients which generate the columns of ${\mathbf{S}}_{\mathit{\text{mn}}}$ as linear combinations of the orthonormal column vectors in such a way that the moment matching conditions are satisfied. The orthonormal column vectors are restricted to the conditions provided in Lemma 1, which helps to reduce the number of equations in the non-linear system. With such construction simulations matching Kollo skewness become highly efficient; moreover, simulations matching Kollo kurtosis become practically possible, but in high-dimensional systems, the kurtosis matrix contains many elements to match, which greatly increases the difficulty of exact simulation.⁷

Our method differs from Alexander et al. (2022) and Hanke et al. (2017), which solve condition (9) for simulations matching up to Kollo skewness. Hanke et al. (2017) set the undetermined elements in the equation system equal to arbitrary values so the equation system becomes a determined system that is potentially solvable. Their method thus requires a trial-and-error process to find arbitrary values, which can be very slow. Moreover, it can fail completely and even when it works the simulations may have undesirable marginal densities for most real-world applications. Noting these issues, Alexander et al. (2022) propose the Random Orthogonal Matrix simulation method with exact Kollo skewness which uses either a statistical bootstrap from real-world data or parametric multivariate distributions as arbitrary values. They also provide the necessary and sufficient conditions on arbitrary values for Kollo skewness matching to be possible.⁸ Such modifications significantly improve the theoretical and computational properties of the simulated samples, but their algorithms require arbitrary values which often produce marginals with strange characteristics. This problem arises because the choice of distribution used to generate random values has a significant impact on the success of solving for the remaining elements. This problem was discussed in depth by Alexander et al. (2022) who show that reducing the variance of the chosen distribution will ensure a low failure rate. But then, by definitions, the remaining elements will have more extreme values so that the first few columns, viewed as samples from marginal distributions of the first few random variables, could be overly leptokurtic. In this case the sample characteristics may not be ideal for many applications, especially when the target values for skewness or kurtosis are large. Our method bypasses this problem completely. With the use of an orthonormal basis, we remove the need to allocate arbitrary values to elements in ${\mathbf{S}}_{\mathit{\text{mn}}}.$ Instead, we can obtain the orthonormal matrices ${\mathbf{\Theta }}_{m,l_k}^k$ in such a way that ${\mathbf{S}}_{\mathit{\text{mn}}}$ has consistent characteristics, which provides a better quality of data, as will be demonstrated below.

Furthermore, direct simulations matching kurtosis become practically impossible when n is large. For instance, when n = 30 the Alexander et al. algorithm requires the solution of 30 different systems of non-linear equations, the smallest being a system of 3 equations (k = 1, i.e., for the first column) and the largest being a system of 32 equations (k = 30, i.e., the last column). In total one needs to solve ${\sum }_{k=1}^n(k+2)=525$ non-linear equations to target Kollo skewness. But if we were to extend the system to target Kollo kurtosis directly, when n = 30 there would be ${\sum }_{k=1}^n(2k+2)=989$ non-linear equations to solve, which may even be beyond the reach of the high-performance computers in use today. Our method provides a practical solution to this problem because the orthonormal basis greatly reduces the number of equations that need to be solved. When matching Kollo skewness, our algorithm only needs to solve one quadratic equation and one cubic (k = 1) or linear ($k=2,\dots ,n$) equation. In addition, we have an analytical solution in Corollary 1 for solving each column of ${\mathbf{S}}_{\mathit{\text{mn}}}$ once the first column of ${\mathbf{S}}_{\mathit{\text{mn}}}$ is available. By contrast, the algorithm proposed by Alexander et al. requires a numerical solution for one quadratic equation and k + 1 linear equations. When matching both Kollo skewness and Kollo kurtosis, our algorithm is much slower because it involves solving two quadratic equations and k linear equations for each column of ${\mathbf{S}}_{\mathit{\text{mn}}}.$ Also, it requires numerical methods to solve a higher-dimensional, non-linear system. Nevertheless, it is the only algorithm that can target the first four multivariate moments exactly.

We finish with a brief demonstration of the flexibility and efficiency of our moment-matching results. Not all matrices ${\mathbf{\Theta }}_{m,l_1}^1$ admit a solution for the conditions of Theorem 1 or 2, so the speed of the general Algorithms 1 and 2 depend on the time it takes to find the matrices ${\mathbf{\Theta }}_{m,l_k}^k$ for which real-valued vectors ${\mathbf{c}}_k$ that satisfy Theorem 1 or 2 exist. Similarly, for $k=2,\dots ,n$ each ${\mathbf{\Theta }}_{m,l_k}^k$ may need to be regenerated several times until we obtain a real-value ${\mathbf{c}}_k$ that satisfies the conditions of Theorem 1 or 2. Although l_k is at least 2 when matching the first three moments, and at least k + 2 when matching the first four moments, if computational speed becomes an issue we prefer to use a greater value for l_k because it is easier to find a real-value solution for an under-determined system than for an exactly determined one. But for matching the first three moments, Corollary 1 provides a simple set of conditions to check if we set l_k = 2, as well as providing the analytical expressions for the elements of ${\mathbf{c}}_k.$ The computation times provided in Figure 1 use l_k = 2 for matching the first three moments and $l_k=k+4$ for matching the first four moments.

Figure 1. Empirical cumulative distribution functions of simulations matching cryptocurrency data. Reported in the left panels, the empirical distribution functions are based on simulations using the algorithm of Alexander et al. (2022) (KROM) and our algorithms (Algorithm 1 for targeting the first three moments and Algorithm 2 for targeting all first four moments) with 1000 observations on three random variables. The target moments are derived from a sample of cryptocurrency data with 720 observations, from 08:00 4 April to 08:00 4 May 2021, UTC time, with Kollo skewness $\mathit{\tau }=[-1.02,-0.03,0.93]\prime$ and the Kollo kurtosis $\mathbf{K}=[7.79,3.97,2.75;3.97,9.30,2.44;2.75,2.44,13.08].$ For simplicity, we standardize all data before setting the targets for simulation and therefore set ${\mathit{\mu }}_n={\mathbf{0}}_n,$ ${\mathbf{\Sigma }}_n={\mathbf{I}}_n$ and ${\mathbf{Q}}_m={\mathbf{I}}_n,$ so that ${\mathbf{X}}_{\mathit{\text{mn}}}={\mathbf{S}}_{\mathit{\text{mn}}}{\mathbf{\Omega }}_n^{\prime }.$ We also employ the rotation matrix shown in Section 3 for all simulations. The black lines in the CDF plots are the empirical distribution functions of the real data on the corresponding variable, fitted using the Gaussian kernel. We treat these lines as the target marginal densities for the KS tests. The example simulations are a typical case of simulation trials that pass the KS test at the significance level of 10%. For simulations using our algorithms, the orthonormal basis is generated using the bootstrapped sample. For simulations using the Alexander et al. algorithm, the arbitrary values in each column are also randomly drawn from the real data and we set the parameter ${\sigma }^2=0.8$ or 0.6. In the right panel, the key characteristics table reports the KS statistics and p-values for the marginals depicted on the left. The table also reports the mean and standard deviation (in brackets) of the average computation time (ACT, in second/simulation), which is computed using the total computational time over 100 trials and depicted in the computational time figure below the table; the red line depicts the computation time of KROM ${}_{0.8}$ and is set against the left axis; the blue and green lines depicts the computation time of KROM ${}_{0.6}$ and our algorithm matching the first three moments and are set against the right axis. The machine used is an Intel Core i5-8250U CPU, 3.41 gigahertz, with 8 gigabyte RAM, running Python under Windows.

Figure 1 presents the empirical comparison of the cumulative density function (CDF) generated by our algorithm for matching the Kollo skewness only, or both the Kollo skewness and kurtosis, against that generated by the algorithm of Alexander et al. (2022).⁹ The Kollo skewness and kurtosis targets are derived empirically from the cryptocurrency hourly returns data. These moments have moderately large elements everywhere, which is suitable to illustrate the two problems discussed above. The figure also reports the results of the Kolmogorov–Smirnov (KS) test which quantifies the distance between the simulated and targeted samples, and the computation time for simulations matching the Kollo skewness only, since we have no benchmark comparison for matching Kollo kurtosis.

First, the results replicate Alexander et al.’s finding that a more constrained selection of arbitrary values (i.e., ${\sigma }^2=0.6,$ in blue, relative to ${\sigma }^2=0.8,$ in red) generates much more leptokurtic samples but is much faster for sample generation. So users must choose between a better-fitted sample or a faster speed of simulation. Note that the example distribution functions reported in Figure 1 are selected from “qualified” trials, i.e., the trials that pass the KS test at the significance level of 10%. Yet we can see the distribution functions generated using Alexander et al.’s algorithm with both ${\sigma }^2=0.6$ and ${\sigma }^2=0.8$ still appear leptokurtic. By comparison, samples generated by our algorithms for matching Kollo skewness and kurtosis do not suffer from the leptokurtic problem. In addition, our algorithm for matching Kollo skewness can be a hundred times faster when targeting suitably distributed samples: it takes 0.08 s per simulation for our algorithm versus 8.04 s per simulation for Alexander et al.’s algorithm with ${\sigma }^2=0.8.$ The computation time for matching all four moments is not reported because it is always much slower than the other algorithms. This is in line with expectations because our algorithm requires at least k + 2 degrees of freedom when solving the kth column of ${\mathbf{S}}_{\mathit{\text{mn}}}.$

Finally, we emphasize that the more multivariate moments that are matched exactly, the more information that is preserved by ROM simulation. Our research thus provides an interesting avenue for future research on matching Kollo kurtosis as well as Kollo skewness—that is to look for an analytical solution for solving each column of ${\mathbf{S}}_{\mathit{\text{mn}}}$ once its first column is available. Similar to the impact of Corollary 1 which significantly speeds up our algorithm for matching Kollo skewness, such a solution should further enhance our algorithm for matching all four moments so that such simulations become not just practically possible, as this paper demonstrates, but much more faster to use in practice.

Correction Statement

This article has been corrected with minor changes. These changes do not impact the academic content of the article.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

1 This category of measures include the Malkovich and Afifi (1973) measures, which are developed specifically for normality tests and these measures are extended by Balakrishnan et al. (2007) and Srivastava (1984); the Isogai (1982) measure which is an extension of the Pearson (1895) mode skewness; the Song (2001) kurtosis measure which is based on Rényi entropy.

2 An exception is the Srivastava (1984) measures which are much easier to compute, but they too cannot be applied to moment-matching simulations. This is because they involve a principal component transformation, which is easy to obtain from a given sample, but it is difficult to generate a sample when we only know the Srivastava moments.

3 The Koizumi et al. measures are defined as the sum of Mardia skewness and kurtosis of all sub-vectors of X_n, but they are complex to compute because the number of sub-vectors increases rapidly with the dimensions X_n.

4 To contrast the number of elements: in a 5-dimensional distribution we have 35 and 70 elements in the full skewness and kurtosis tensors, but only 5 and 15 elements in Kollo’s measures.

5 See Ledermann et al. (2011) for discussions on the usual conditions for rotation matrices. Ledermann and Alexander (2012) further investigate the sampling properties of different types of rotation matrices.

6 Alexander et al. (2022) note a similar property for Kollo skewness.

7 In both cases our method may be regarded as a generalization of the original ROM simulation to higher-moment matching, because using only one orthonormal vector to generate each column of ${\mathbf{S}}_{\mathit{\text{mn}}}$ reduces our method to a ROM simulation which only targets the mean vector and covariance matrix.

8 If the constraints are not satisfied, the arbitrary values require re-generation until the equation system becomes solvable with these values.

9 We do not include the algorithm of Hanke et al. (2017) because Alexander et al. (2022) have already demonstrated that the modifications provided by their algorithm produce marginals with superior characteristics and it is very much more computationally efficient. The Hanke et al. method may also suffer from a high failure rate.

Appendix A.

Real-world examples

We illustrate the Kollo skewness and kurtosis and compare them with the popular Mardia measures using a set of three hourly cryptocurrency returns (bitcoin, ether, and litecoin) and the eleven different U.S. stock sector index daily returns, from which we show just three sectors (energy, finance, and real estate) although the Kollo measures are computed from the entire eleven-dimensional distribution. The returns are presented in Figure A1.

Figure A1. Time series of financial returns. (a) Exhibits the hourly returns on three cryptocurrencies: bitcoin (BTC), ether (ETH), and litecoin (LTC) from 1 January 2017 to 31 December 2021; and (b) exhibits the daily returns for three U.S. stock sector indices, viz. Energy, Finance, and Real Estate, from 16 October 2001 to 31 December 2021.

Given that each series is clearly far from i.i.d., we used a fixed-size rolling window to estimate the sample Kollo skewness and kurtosis. The sample size is 720 (∼1 month) for the hourly cryptocurrency data and 250 (∼ one year) for the daily stock index data. The relevant elements of the Kollo skewness vector and kurtosis matrix are depicted as time series in Figure A2.

Figure A2. Kollo skewness and kurtosis in financial data in a rolling-window framework. We plot all elements in the Kollo skewness vector and only plot the unique elements in the Kollo kurtosis matrix. (a) Exhibits the Kollo skewness and kurtosis using a rolling window of 720 hourly returns for cryptocurrency prices of bitcoin (BTC), ether (ETH), and litecoin (LTC) from 1 January 2017 to 31 December 2021. (b) Exhibits the Kollo skewness and kurtosis using a rolling window of 250 daily returns for three U.S. stock sector indices, viz. Energy, Finance, and Real Estate, from 16 October 2001 to 31 December 2021. The grey dotted line depicts the rolling Mardia skewness and kurtosis, using the right-hand scale.

The spikes in skewness and kurtosis clearly represent significant events in each sector. For example, during the financial crisis of 2008, we see a spike in both τ₇ and ${\kappa }_{7,7}$ following the collapse of Lehman Brothers in September 2008. The extremely low values of τ₁ and the spike in ${\kappa }_{1,1}$ during the first half of 2018 were caused by an upward surge in oil prices. Similarly, the upward spike in τ₇ and in most elements in Kollo kurtosis on 16 March 2020 resulted from extreme negative returns across all financial markets following the outbreak of Covid-19.

Almost all elements of the Kollo skewness vector lie in the range $[-10,10]$ with the most sampling variability in τ₁ and τ₃. The Kollo kurtosis ranges from 0 to 200. The spikes in Kollo skewness and Kollo kurtosis occur simultaneously. In particular, the global crash in risky assets after the outbreak of Covid in January 2020 is evident from the dip in Kollo skewness down to almost −10 in both data sets. The Kollo kurtosis takes a narrower range in the stock sector data than for the cryptocurrency data.

Compared with the Mardia measures, we see that the Kollo measures are more sensitive to the extremes of the distribution. Unlike the Mardia skewness, the Kollo measure also captures the direction of asymmetry in the marginals. The Kollo measures also capture more information from co-moments and not only from the marginals. They can exhibit substantial variability even when the Mardia measures remain relatively stable over time. Furthermore, the Kollo measures can indicate which variable contributes most to an extreme value for skewness and kurtosis but Mardia measures only capture the overall information, in a scalar value.

Appendix B.

Mathematical proofs

B.1 Proof of Lemma 1

Write ${\mathbf{\Theta }}_{m,l_k}^k=\left[{\mathit{\theta }}_1^k,\dots ,{\mathit{\theta }}_{l_k}^k\right],$ and ${\mathbf{c}}_k=\left[c_1,\dots ,c_{l_k}\right]\prime$ where at least one element is non-zero. Now if the column vectors in ${\mathbf{\Theta }}_{m,l_k}^k$ are also orthogonal to the vectors ${\mathbf{s}}_1^c,\dots ,{\mathbf{s}}_{k-1}^c$ we have: $$\hspace{1em}\mathbf{1}{\prime }_m({\mathit{s}}_h^c°{\mathit{\theta }}_j^k)=0\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}j=1,\dots ,l_k\hspace{1em}\text{and}\hspace{1em}\forall h<k.$$

Then, since ${\mathbf{\Theta }}_{m,l_k}^k$ is orthogonal to ${\mathbf{1}}_m$ we have ${\mathbf{1}}_m^{\prime }{\mathit{\theta }}_i^k=0$ and hence: $${\mathbf{1}}_m^{\prime }{\mathit{s}}_k^c={\mathbf{1}}_m^{\prime }\left({\mathbf{\Theta }}_{m,l_k}^k {\mathit{c}}_k\right)={\mathbf{1}}_m^{\prime }\sum _{i=1}^{l_k}{c}_i{\mathit{\theta }}_i^k=\sum _{i=1}^{l_k}{c}_i{\mathbf{1}}_m^{\prime }{\mathit{\theta }}_i^k=0,$$ and for $j<k$ $$\hspace{1em}{\mathbf{1}}_m^{\prime }({\mathit{s}}_j^c°{\mathit{s}}_k^c)={\mathbf{1}}_m^{\prime }[{\mathit{s}}_j^c°\left({\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k\right)]={\mathbf{1}}_m^{\prime }({\mathit{s}}_j^c°\sum _{i=1}^{l_k}{c}_i{\mathit{\theta }}_i^k)=\sum _{i=1}^{l_k}{c}_i{\mathbf{1}}_m^{\prime }({\mathit{s}}_j^c°{\mathit{\theta }}_i^k)=0.$$

B.2 Proof of Theorem 1

When k = 1 the conditions for ${\mathbf{s}}_1^c$ are set out in Equation (11)(11) $$\sum _{i=1}^m{s}_{i1}=0\hspace{0.17em},\hspace{1em}\sum _{i=1}^m{s}_{i1}{}^2=m\hspace{0.17em},\hspace{1em}\sum _{i=1}^m{s}_{i1}{}^3=m{\tilde{\tau }}_1.$$(11) and by Lemma 1 we have ${\mathbf{s}}_1^c={\mathbf{\Theta }}_{m,l_k}^1{\mathbf{c}}_1$ so we only need to solve the second two equations in (11)(11) $$\sum _{i=1}^m{s}_{i1}=0\hspace{0.17em},\hspace{1em}\sum _{i=1}^m{s}_{i1}{}^2=m\hspace{0.17em},\hspace{1em}\sum _{i=1}^m{s}_{i1}{}^3=m{\tilde{\tau }}_1.$$(11) . The quadratic equation can be written as: $$m={\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c)={{\mathit{s}}_1^c}^{\prime }{\mathit{s}}_1^c={\left({\mathbf{\Theta }}_{m,l_k}^1{\mathit{c}}_1\right)}^{\prime }\left({\mathbf{\Theta }}_{m,l_k}^1{\mathit{c}}_1\right)={\mathit{c}}_1^{\prime }{{\mathbf{\Theta }}_{m,l_k}^1}^{\prime }{\mathbf{\Theta }}_{m,l_k}^1{\mathit{c}}_1={\mathit{c}}_1^{\prime }{\mathit{c}}_1,$$ because ${{\mathbf{\Theta }}_{m,l_k}^1}^{\prime }{\mathbf{\Theta }}_{m,l_k}^1={\mathbf{I}}_n.$ And the cubic equation becomes: $$m{\tilde{\tau }}_1={\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_1^c)={\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)],$$ which proves (13a)(13a) $$\text{for}\hspace{1em}k=1:{\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_k}^1 {\mathit{c}}_1)]=m{\tilde{\tau }}_1\hspace{0.17em},$$(13a) . When $k=2,\dots ,n$ we need to solve a quadratic Equation (9c)(9c) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_k^c°{\mathit{s}}_k^c)=m\hspace{0.17em},$$(9c) and a linear Equation (9d)(9d) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_k^c)=m{\tilde{\tau }}_k\hspace{0.17em},$$(9d) . The quadratic equation may be written as: $$m={\mathbf{1}}_m^{\prime }({\mathit{s}}_k^c°{\mathit{s}}_k^c)={{\mathit{s}}_k^c}^{\prime }{\mathit{s}}_k^c={\mathit{c}}_k^{\prime }{{\mathbf{\Theta }}_{m,l_k}^k}^{\prime }{\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k={\mathit{c}}_k^{\prime }{\mathit{c}}_k\hspace{0.17em},$$ as before. The linear equation becomes: $$\begin{array}{c}m{\tilde{\tau }}_{}={\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_k^c)={({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{s}}_k^c={({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }\left({\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k\right)\\ =\left({({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathbf{\Theta }}_{m,l_k}^k\right){\mathit{c}}_k=[{({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{\theta }}_1^k,\dots ,{({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{\theta }}_{l_k}^k]{\mathit{c}}_k\\ ={\mathbf{1}}_m^{\prime } [{\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{\theta }}_1^k,\dots ,{\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{\theta }}_{l_k}^k] {\mathit{c}}_k={\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k\hspace{0.17em},\end{array}$$ which proves (13b)(13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b) .

B.3 Proof of Corollary 1

Equation (12)(12) $${\mathit{c}}_k^{\prime }{\mathit{c}}_k=m\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for}\hspace{1em}k=1,\dots ,n;$$(12) defines a 2-dimensional sphere and Equation (13b)(13b) $$\text{for}\hspace{1em}k=2,\dots ,n:{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\tau }}_k\hspace{0.17em},$$(13b) defines a 2-dimensional plane. Hence the question of whether the system of equations has real roots can be converted into the question of whether the sphere has intersections with the plane. We know that a necessary and sufficient condition for this is that the shortest distance between the centre of the sphere and the plane should be no more than the radius of the sphere, i.e., $$\frac{\left|m{\tilde{\tau }}_k\right|}{\sqrt{{\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,2}^k{{\mathit{B}}_{m,2}^k}^{\prime }{\mathbf{1}}_m}}\le \sqrt{m}\hspace{0.17em},$$ which proves (14)(14) $${\mathbf{1}}_m^{\prime }{\mathit{B}}_{m,2}^k {{\mathit{B}}_{m,2}^k}^{\prime }{\mathbf{1}}_m\ge m{{\tilde{\tau }}_k}^2\hspace{0.17em},$$(14) . The sphere and the linear equation for the plane may be written as follows, $$\begin{array}{c}{\left(c_1^k\right)}^2+{\left(c_2^k\right)}^2=m\\ {({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{\theta }}_1^k{c}_1^k+{({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{\theta }}_2^k{c}_2^k=m{\tilde{\tau }}_k.\end{array}$$

Hence we may set: (B1) $$c_2^k=\frac{m{\tilde{\tau }}_k-{({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{\theta }}_1^k{c}_1^k}{{({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{\theta }}_2^k},$$(B1) and so the sphere equation becomes: $${\left(c_1^k\right)}^2+{\left(\frac{m{\tilde{\tau }}_k-{({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{\theta }}_1^k{c}_1^k}{{({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{\theta }}_2^k}\right)}^2=m.$$

Let ${\tilde{b}}_1^k$ and ${\tilde{b}}_2^k$ denote the first and second elements of ${\mathbf{1}}_m^{\prime }{\mathbf{B}}_{m,2}^k,$ i.e., ${\tilde{b}}_1^k=({\mathbf{s}}_1^c°{\mathbf{s}}_1^c)\prime {\mathit{\theta }}_1^k$ and ${\tilde{b}}_2^k={({\mathbf{s}}_1^c°{\mathbf{s}}_1^c)}^{\prime }{\mathit{\theta }}_2^k.$ Then the above equation becomes: $$\left[{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2\right]{\left(c_1^k\right)}^2-2b_1^{k}m{\tilde{\tau }}_k{c}_1^k+m^2{\tilde{\tau }}_k^{ 2}-m{\left({\tilde{b}}_2^k\right)}^2=0.$$

Using the standard formula for the roots of this equation, the solutions are: $$\begin{array}{c}{c}_1^k=\frac{2{\tilde{b}}_1^{k}m{\tilde{\tau }}_k\pm \sqrt{4{\left({\tilde{b}}_1^k\right)}^2{m}^2{\tilde{\tau }}_k{}^{ 2}-4\left[{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2\right]\left(m^2{\tilde{\tau }}_k{}^{ 2}-m{\left({\tilde{b}}_2^k\right)}^2\right)}}{2\left[{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2\right]}\\ =\frac{{\tilde{b}}_1^{k}m{\tilde{\tau }}_k\pm \sqrt{\left[{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2\right]{\left({\tilde{b}}_2^k\right)}^{2}m-{\left({\tilde{b}}_2^k\right)}^2{m}^2{\tilde{\tau }}_k{}^2}}{{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2}\\ =\frac{{\tilde{b}}_1^{k}m{\tilde{\tau }}_k\pm {\tilde{b}}_2^k\sqrt{\left[{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2\right]m-m^2{\tilde{\tau }}_k{}^2}}{{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2}.\end{array}$$

Finally, we substitute these values for $c_1^k$ into (B1)(B1) $$c_2^k=\frac{m{\tilde{\tau }}_k-{({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{\theta }}_1^k{c}_1^k}{{({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{\theta }}_2^k},$$(B1) to obtain: $$c_2^k=\frac{m{\tilde{\tau }}_k-{({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{\theta }}_1^k{c}_1^k}{{({\mathit{s}}_1^c°{\mathit{s}}_1^c)}^{\prime }{\mathit{\theta }}_2^k}=\frac{{\tilde{b}}_2^{k}m{\tilde{\tau }}_k\mp {\tilde{b}}_1^k\sqrt{[{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2]m-m^2{{\tilde{\tau }}_k}^2}}{{\left({\tilde{b}}_1^k\right)}^2+{\left({\tilde{b}}_2^k\right)}^2}.$$

B.4 Proof of Theorem 2

To match Kollo kurtosis each column of ${\mathbf{S}}_{\mathit{\text{mn}}}$ needs to satisfy Equation (9e)(9e) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c°{\mathit{s}}_k^c)=m{\tilde{\kappa }}_{\mathit{\text{ik}}},\forall i\le k,$$(9e) in addition to satisfying the conditions (9a)—(9d), thus (i)—(iii) has been proven in Theorem 1 and the following focuses on (iv). When k = 1, ${\mathbf{s}}_1^c$ needs to satisfy a quartic Equation (16)(16) $$\sum _{i=1}^m{s}_{i1}{}^4=m{\tilde{\kappa }}_{11}\hspace{0.17em},$$(16) which is easily rewritten as (17a)(17a) $${\mathbf{1}}_m^{\prime }[({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)°({\mathbf{\Theta }}_{m,l_1}^1 {\mathit{c}}_1)]=m{\tilde{\kappa }}_{11}\hspace{0.17em},$$(17a) , viz. $$m{\tilde{\kappa }}_{11}={\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_1^c)={\mathbf{1}}_m^{\prime }[\left({\mathbf{\Theta }}_{m,l_1}^1{\mathit{c}}_1\right)°\left({\mathbf{\Theta }}_{m,l_1}^1{\mathit{c}}_1\right)°\left({\mathbf{\Theta }}_{m,l_1}^1{\mathit{c}}_1\right)°\left({\mathbf{\Theta }}_{m,l_1}^1{\mathit{c}}_1\right)].$$

When $k=2,\dots ,n,$ ${\mathbf{s}}_k^c$ should satisfy condition (9e)(9e) $${\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c°{\mathit{s}}_k^c)=m{\tilde{\kappa }}_{\mathit{\text{ik}}},\forall i\le k,$$(9e) , which contains a quadratic equation when i = k and k − 1 linear equations for $i=1,\dots .,k-1.$ First, we focus on the quadratic equation: $$\begin{array}{c}m{\tilde{\kappa }}_{\mathit{\text{kk}}}={\mathbf{1}}_m^{\prime }\left({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_k^c°{\mathit{s}}_k^c\right)=\left({\mathit{s}}_1^c°{\mathit{s}}_k^c\right)\prime \left({\mathit{s}}_1^c°{\mathit{s}}_k^c\right)\hfill \\ =\left({\mathit{s}}_1^c°\left({\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k\right)\right)\prime \left({\mathit{s}}_1^c°\left({\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k\right)\right)\hfill \\ ={\mathit{c}}_k^{\prime }\left[{\mathit{s}}_1^c°{\mathit{\theta }}_1,\dots ,{\mathit{s}}_1^c°{\mathit{\theta }}_{l_k}\right]\prime \left[{\mathit{s}}_1^c°{\mathit{\theta }}_1,\dots ,{\mathit{s}}_1^c°{\mathit{\theta }}_{l_k}\right]{\mathit{c}}_k\hfill \\ ={\mathit{c}}_k^{\prime }{\mathit{D}}_{m,l_k}^k\prime {\mathit{D}}_{m,l_k}^k{\mathit{c}}_k\hspace{0.17em},\hfill \end{array}$$ which proves the condition (17b)(17b) $${\mathit{c}}_k^{\prime }{{\mathit{D}}_{m,l_k}^k}^{\prime }{\mathit{D}}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\kappa }}_{\mathit{\text{kk}}}\hspace{0.17em},$$(17b) . Next, consider these k − 1 linear equations, for each i < k, $$m{\tilde{\kappa }}_{\mathit{\text{ik}}}={\mathbf{1}}_m^{\prime }({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c°{\mathit{s}}_k^c)={({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c)}^{\prime }{\mathit{s}}_k^c={({\mathit{s}}_1^c°{\mathit{s}}_1^c°{\mathit{s}}_i^c)}^{\prime }\left({\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k\right).$$

Then the Equation (17c)(17c) $$\text{and}\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}{{\mathit{H}}_{m,k-1}^k}^{\prime }{\mathbf{\Theta }}_{m,l_k}^k{\mathit{c}}_k=m{\tilde{\mathit{\kappa }}}_{k-1},$$(17c) is simply rewritten in a matrix form.