Inducing a desired value of correlation between two point-scale variables: a two-step procedure using copulas

Focusing on point-scale random variables, i.e. variables whose support consists of the first m positive integers, we discuss how to build a joint distribution with pre-specified marginal distributions and Pearson’s correlation ρ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document}. After recalling how the desired value ρ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} is not free to vary between -1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1$$\end{document} and +1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+1$$\end{document}, but generally ranges a narrower interval, whose bounds depend on the two marginal distributions, we devise a procedure that first identifies a class of joint distributions, based on a parametric family of copulas, having the desired margins, and then adjusts the copula parameter in order to match the desired correlation. The proposed methodology addresses a need which often arises when assessing the performance and robustness of some new statistical technique, i.e. trying to build a huge number of replicates of a given dataset, which satisfy—on average—some of its features (for example, the empirical marginal distributions and the pairwise linear correlations). The proposal shows several advantages, such as—among others—allowing for dependence structures other than the Gaussian and being able to accommodate the copula parameter up to an assigned level of precision for ρ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} with a very small computational cost. Based on this procedure, we also suggest a two-step estimation technique for copula-based bivariate discrete distributions, which can be used as an alternative to full and two-step maximum likelihood estimation. Numerical illustration and empirical evidence are provided through some examples and a Monte Carlo simulation study, involving the CUB distribution and three different copulas; an application to real data is also discussed.


Introduction
Datasets arising in the social sciences often contain ordinal variables. Sometimes they are genuine ordered assessments (judgements, preferences, degree of liking of a product or adhesion to a sentence, etc.), whereas in other circumstances they are discretized or categorized for convenience (age of people in classes, education achievement, levels of blood pressure, etc.) The former situation often arises when a survey is administered to a group of people being studied, e.g. questionnaires submitted by a company to their customers with the aim of assessing their level of satisfaction towards a product or service the company has provided. Respondents choose a qualitative assessment on a graduated sequence of verbal definitions (for instance, "extremely dissatisfied", "very dissatisfied" , … "very satisfied", "extremely satisfied"), also known as "Likert scale", which can be coded as integers numbers ( 1, 2 … , m ) just for convenience: this amounts to assuming that the categories are evenly spaced (Iannario and Piccolo 2012). There are several statistical models and techniques that can be employed for handling multivariate ordinal data without trying to quantify their ordered categories. (The review by Liu and Agresti (2005) and the later textbook of Agresti (2010) give a thorough treatment.) Among them, correlation models and association models both study departures from independence in contingency tables and involve the assignment of scores to the categories of the row and column variables in order to maximize the relevant measure of relationship (the correlation coefficient in the correlation models or the measure of intrinsic association in association models, see Faust and Wasserman 1993). We also remind nonlinear principal component analysis (NLPCA), which is a special case of a multivariate reduction technique named homogeneity analysis and which can be usefully applied in customer satisfaction surveys (Ferrari and Manzi 2010) for mapping the observed ordinal variables into a one-dimensional (or, more generally, lower-dimensional) quantitative variable. Neither the weights of the original variables nor the differences between their categories (this is the distinction with standard principal component analysis) are assumed a priori.
However, substituting the ordered categories with the corresponding integer numbers, though representing just an arbitrary assumption, is still quite a common and accepted practice, which leads to further multivariate statistical analyses handling them as (correlated) discrete variables (Norman 2010;Carifio and Perla 2008). Now, one may be interested in building and simulating a multivariate random vector whose univariate components are point-scale variables. In fact, describing a real phenomenon by creating mirror images and imperfect proxies of the (partially) unknown underlying population in a repeated manner allows researchers to study the performance of their statistical methods through simulated data replicates that mimic the real data characteristics of interest in any given setting (Demirtas and Yavuz 2015;Demirtas and Vardar-Acar 2017). This is often necessary since exact analytic results are seldom available for finite sample sizes, and thus, simulation is required to assess the reliability, validity, and plausibility 1 3 Inducing a desired value of correlation between two point-scale… of inferential techniques and to evaluate to which extent they are robust to deviations from statistical assumptions. However, rather than completely detailing a joint distribution for modelling the phenomena under study, it is often more convenient and realistic to specify only the marginal distributions and pairwise correlations, which are very easy to interpret and whose sample analogues can be easily computed on the dataset at hand one wants to "reproduce". In Lee (1997), some methods are described for generating random vectors of categorical or ordinal variables with specified marginal distributions and degrees of association between variables. For ordinal variables, a common index for measuring association is Goodman and Kruskal's coefficient (Kruskal and Goodman 1954;Ruiz and Hüllermeier 2012), ranging between −1 and +1 , with zero corresponding to independence. A first proposal is based upon using convex combination of joint distribution with extremal values of (extremal tables); a second one relies on threshold arguments and involves Archimedean copulas. In this paper, we suggest a sort of modification of this latter method to correlated point-scale variables.
In the following, we will limit our analysis to the bivariate case, which is by far easier to deal with, but whose results, with some caution, can be generalized to the multivariate context. We consider two point-scale random variables (rvs), X 1 and X 2 , defined over the support spaces X 1 = 1, 2, … , m 1 and X 2 = 1, 2, … , m 2 , respectively, with probability mass functions p 1 (i) = p i⋅ = P(X 1 = i), i = 1, … , m 1 , and p 2 (j) = p ⋅j = P(X 2 = j), j = 1, … , m 2 . We want to determine some bivariate probability mass function p ij = p(i, j) = P(X 1 = i, X 2 = j), i = 1, … , m 1 ;j = 1, … , m 2 such that its margins are p 1 and p 2 and the correlation X 1 ,X 2 is equal to an assigned . In order to give an answer to this question, we have first to recall two properties of Pearson's correlation, which apply to both the continuous and, to even a larger extent, the discrete case; this is the topic of Sect. 2. In Sect. 3, we first state the problem of finding a joint probability function with assigned margins and correlation in general terms; then, we focus on a particular class of joint distributions, recalling how to build copula-based bivariate discrete distributions; finally, we describe the proposed procedure for inducing a desired value of correlation between two point-scale variables. Section 4 illustrates an application to CUB distributions. Section 5 recalls inferential procedures for dependent rvs and based on the algorithm of Sect. 3 devises a sort of moment method for estimating the dependence parameter of the copula-based bivariate distribution. Section 6 describes a Monte Carlo simulation study whose aim is to comparatively assess the statistical performances of the new inferential method and the existing ones based on maximum likelihood. Section 7 provides an application to a real data set. In the concluding section, some final remarks are provided.

Attainable correlations between two random variables
Pearson's linear correlation is by far the most popular measure of correlation between two quantitative variables.
It is often employed to measure correlation also between Likert scale variables, which arises some criticism.

3
Despite the many useful properties it enjoys, it also reveals some disadvantages. A first drawback of Pearson's correlation is that given two marginal cumulative distribution functions (cdfs) F 1 and F 2 and a correlation value ∈ [−1, +1] , it is not always possible to construct a joint distribution F with margins F 1 and F 2 , whose correlation is equal to the assigned . This is an issue that is often underrated if not neglected by researchers (Leonov and Qaqish 2020). We can state this result (often reported as "attainable correlations", see McNeil et al. 2005, pp. 204-205) in the following way. Let (X 1 , X 2 ) be a random vector with marginal cdfs F 1 and F 2 and an unspecified joint cdf; assume also that Var(X 1 ) > 0 and Var(X 2 ) > 0 . The following statements hold: 1. The attainable correlations form a closed interval [ min , max ] with min < 0 < max . 2. The minimum correlation = min is attained if and only if X 1 and X 2 are countermonotonic. The maximum correlation = max is attained if and only if X 1 and X 2 are comonotonic. 3. min = −1 if and only if X 1 and −X 2 are of the same type, and max = 1 if and only if X 1 and X 2 are of the same type.
For point-scale rvs X 1 and X 2 , it is then clear that the maximum correlation is +1 if and only if they are identically distributed: In the general case, the values min and max can be computed by building the cograduation and countergraduation tables (see Salvemini 1939 andBarbiero 2012 for an example of calculation).
A second result about Pearson's correlation can be resumed as follows: given two margins F 1 and F 2 and a feasible linear correlation (i.e. a value falling within the interval [ min , max ] ), the joint distribution F having margins F 1 and F 2 and correlation is not unique. In other terms, the marginal distributions and pairwise correlation of a bivariate rv do not univocally determine its joint distribution. Even if this second fallacy may represent a limit from one side, on the other side represents a form of flexibility, since it means that given two pointscale distributions and a consistent value of , there are generally several (possibly, infinite) different ways to join them into a bivariate distribution with that value of correlation, as we will see in the next two sections.

A procedure for inducing a desired value of correlation between two point-scale random variables with assigned marginal distributions
We will now state the problem object of this work in general terms; then, resorting to copulas, we will reformulate it in a more specific context. Somehow, we will split the original problem into two sequential sub-problems: (i) finding a 1 3 Inducing a desired value of correlation between two point-scale… family of joint distributions with the assigned margins, (ii) finding within this family a distribution with the desired value of correlation.

Statement of the problem
The problem of finding a bivariate point-scale distribution with assigned marginal distributions and Pearson's correlation can be laid out as follows. We have to find the m 1 × m 2 probabilities p ij , 0 ≤ p ij ≤ 1 , defining the joint pmf of the rv (X 1 , X 2 ) , which satisfy the following system of equalities: The first two equalities correspond to the request of matching the assigned marginal distributions and the latter to the assigned correlation. The total number of equality constraints is m 1 + m 2 ( m 1 + m 2 − 1 actual constraints on the two margins, plus one on Pearson correlation).
If, for example, m 1 = m 2 = 2 (i.e. if X 1 and X 2 are shifted Bernoulli rvs), then we have a system of 4 equations in 4 variables, which yields a unique solution for the p ij . In this case, in fact, one can easily prove that if the assigned falls within the bounds min and max , equal to the probabilities satisfying the system (1) are For higher values of m 1 (and m 2 ) the solution is not unique; generally, there are infinite solutions (i.e. bivariate distributions) satisfying system (1)-given that the bivariate correlation bounds are respected (i.e. min ≤ ≤ max ). We illustrate it through the following example.
Instead of considering the whole set of feasible joint probabilities p ij obtained by solving system (1), one can restrict the analysis to a particular subset of these solutions satisfying the first two constraints on margins. For this aim, we will now recall the concept of copula and copula-based bivariate discrete distributions.

Generating bivariate discrete distributions having the pre-specified margins
Using copulas represents a straightforward solution for easily constructing a multivariate distribution respecting the assigned margins. A d-dimensional copula C is a joint cdf in [0, 1] d with standard uniform margins U j , j = 1 … , d: The importance of copulas in studying multivariate cdfs is summarized by the Sklar's theorem (McNeil et al. 2005), whose version for d = 2 states that if F 1 and F 2 are the cdfs of the rvs X 1 and X 2 , the function defines a valid joint cdf, whose margins are exactly F 1 and F 2 . This result keeps holding if X 1 and X 2 are point-scale rvs; in this case, the joint pmf can be derived from (2) as: for i = 1, … , m 1 ;j = 1, … , m 2 . It is worth noting that given a joint cdf F with margins F 1 and F 2 , Sklar's theorem also states that there exists a copula C such that F can be written as in (2). This copula is unique if F 1 and F 2 are continuous; on the contrary, uniqueness is not guaranteed if they are discrete.
There exists a multitude of copulas, which-as happens for joint cdfs-usually depend on some parameter ; we will now review three well-known parametric copula families.

The Gauss copula
The d-variate Gauss copula is the copula that can be extracted from a d-variate normal vector Y Y Y with mean vector and covariance matrix and is exactly the same as the copula of X X X ∼ N d (0 0 0, P) , where P is the correlation matrix of Y Y Y . In two dimensions, it can be expressed, for Ga ≠ ±1 , as: ds 1 ds 2 .

The Frank copula
The one-parameter bivariate Frank copula is defined as with ≠ 0 . For → 0 , we have that the Frank copula reduces to the independence copula; for → +∞ , it tends to the comonotonicity copula; for → −∞ , it tends to the countermonotonicity copula.

The Plackett copula
The one-parameter bivariate Plackett copula is defined as with ∈ (0, +∞) ⧵ {1} . When → 1 , it reduces to the independence copula, whereas for → 0 it tends to the countermonotonicity copula and for → ∞ to the comonotonicity copula.

An algorithm for inducing any feasible value of correlation within a parametric copula-based family of distributions
In order to induce any feasible value of correlation between the two discrete margins of the distribution (2), we have further to impose that the copula C(⋅; ) is able to encompass the entire range of dependence, from perfect negative dependence (which leads to the linear correlation min ) to perfect positive dependence ( max ). Copulas enjoying this property are named "comprehensive"; the three copulas recalled in the previous section are all comprehensive.
Once the marginal distributions of X 1 and X 2 are assigned, and the parametric copula C(⋅; ) has been selected, their correlation coefficient X 1 ,X 2 will depend only on the copula parameter ∈ [ min , max ] ; this relationship may be written in an analytical or numerical form, say X 1 ,X 2 = g( |F 1 , F 2 ) . Since the function g is not usually analytically invertible, inducing a desired value of correlation between two point-scale variables, falling in [ min , max ] , by setting an appropriate value of , is a task that can be generally done only numerically, by finding the (unique) root of the equation g( ) − = 0 . If g( ) is a monotone increasing function of the copula parameter, it can be implemented by resorting to the following iterative procedure (see a similar proposal for the Gauss copula in Ferrari and Barbiero 2012; Barbiero and Ferrari 2017; and an early extension to other copulas in Barbiero 2018): 1. Set (0) = (with being the value of for which the copula C reduces to the independence copula); (0) = 0. 2. Set t ← 1 and = (t) , with (t) some value strictly greater (smaller) The iterative process at the basis of the algorithm is quite clear (see Fig. 1): one starts from two points in the ( , ) Cartesian diagram: A = ( (0) , (0) ) and B = ( (1) , (1) ) , where (1) can be chosen arbitrarily, respecting the unique condition that the resulting (1) has the same sign as the target . From these two points, one derives the next value of , (2) (corresponding to the abscissa of point C), by linear interpolation, considering the slope m (2) , associated with the line passing through them, respecting the lower or upper bounds of (this is why the min and max operators appear in the recursive formulas for (t) of step 6); the procedure then continues, computing the actual value (2) (ordinate of point D), and then iteratively updating (t) (and computing (t) ) by taking into account just the last two points for determining the updated m (t) .
The above heuristic algorithm makes sense if g is a monotone increasing function, which is often the case: for the Gauss, Frank, and Plackett copulas, the linear correlation is an increasing function of the dependence parameter , keeping fixed the two marginal distributions. In fact, let us recall that we say that the joint cdf F(x 1 , x 2 ; ) , with fixed margins F 1 and F 2 , is "increasing in concordance" as increases if, for any 2 > 1 , Then, it follows (see, for example, Scarsini and Shaked 1996) Since the Gauss copula, the Frank copula, and the Plackett copula are all increasing in concordance with respect to their parameter, i.e. (Joe 2014) and the same holds for the joint cdf F(x 1 , x 2 ; ) = C(F 1 (x 1 ), F 2 (x 2 ); ) , then we can claim that X 1 ,X 2 is increasing in .
A very particular case is represented by the Gauss copula. A known theoretical result, which goes under the name of Lancaster theorem (Lancaster 1957, p.290), and later reported for example in Cario and Nelson (1997), allows us to claim that the correlation between the discrete rvs X 1 and X 2 has the same sign of and in absolute value it is not greater than the Gauss copula correlation: sgn( X 1 ,X 2 ) = sgn( Ga ) and | X 1 ,X 2 | ≤ | Ga | . Therefore, a reasonable value for the starting value (1) ∶= (1) Ga is the value of the target correlation itself.
The advantage of the proposed algorithm stands in the four following (connected) features: (i) the flexibility in the choice of the underlying copula, which can be different from the Gaussian and is just required to span the entire dependence range, (ii) the capacity of finding the appropriate value of without making use of any sample from the two marginal distributions, thus avoiding introducing sampling errors, (iii) the possibility of controlling a priori the error (absolute difference between target and actual value of X 1 ,X 2 )-setting equal to 10 −7 generally allows to recover in a few steps, and (iv) the absence of inner potentially time-consuming optimization or root finding routines.
Existing procedures for solving the same (or a similar) problem are available in the literature, but do not enjoy all the features above mentioned. For example, the proposal by Demirtas (2006) is based on simulating binary data whose marginals are derived collapsing the prespecified marginals of ordinal variables. The correlation matrix of the binary variables is obtained by an iterative process in order to match the target correlation matrix for ordinal data, which requires the generation of a "huge" bivariate sample of binary data.
Other proposals by Madsen and Dalthorp (2007), Ferrari and Barbiero (2012), or Xiao (2017) for simulating multivariate ordinal/discrete variables with assigned margins and correlations exclusively address the dependence structure induced by the Gauss copula. Lee and Kaplan (2018) proposed two procedures based on the principles of maximum entropy and minimum cross-entropy to simulate multivariate ordinal variables with assigned values of marginal skewness and kurtosis; they rely on the multivariate normal distribution as a latent variable. Foldnes and Olsson (2016) proposed a simulation technique for nonnormal data with pre-specified skewness, kurtosis, Inducing a desired value of correlation between two point-scale… and covariance matrix, by using linear combinations of independent generator (IG) variables; its most important feature is that the resulting copula is not Gaussian. In Nelsen (1987), using convex linear combinations of the pmfs for the discrete Fréchet boundary distributions (i.e. those corresponding to comotonicity and countermonotonicity) and the pmf for independent rvs, the author constructs bivariate pmfs for dependent discrete rvs with arbitrary marginals and any correlation between the theoretical minimum and maximum. A similar rationale has been later used by Demirtas and Vardar-Acar (2017) for devising an algorithm for inducing any desired Pearson or Spearman correlation to independent bivariate data whose marginals can be of any distributional type and nature. The algorithm we proposed, though limited to point-scale rvs, is much more flexible as it allows a much broader choice of dependence structures, whereas the latter two procedures employ a convex combination of bivariate comonotonicity and countermonotonicity copulas. (An analogous way was followed by Lee (1997) for the first method.) Obviously, in the simplest case of a bivariate shifted Bernoulli ( m 1 = m 2 = 2 ), the proposed algorithm recovers (numerically) the same unique bivariate distribution yielded (analytically) by system (1), whatever copula is selected (provided it spans the entire dependence spectrum, i.e. it is comprehensive). For the case of dependent Bernoulli rvs, see also the example presented in McNeil et al. 2005, p.188.
We remark that even if we mentioned three well-known comprehensive copulas (Gauss, Frank, and Plackett), which are all exchangeable and radially symmetric (see, for example, McNeil et al. (2005), chapter 5), exchangeability and radial symmetry are not necessary conditions for the algorithm to work; thus, the proposed procedure is able to deal with asymmetrical dependence, which often occurs in many fields, especially in finance. The comprehensive property is instead required if one wants to span the entire range of feasible linear correlations between the two point-scale rvs. If one uses the Gumbel or the Clayton copula (both belonging to the broad class of Archimedean copulas) to induce correlation between the rvs, since these two copulas can only model positive dependence through their scalar parameter , it descends that only positive (or at most null) values of linear correlation can be induced and then assigned. A useful reference is Table 4.1 in Nelsen (2006), where some important one-parameter families of Archimedean copulas are listed along with some special and limiting cases; from there, it is thus possible to distinguish and select comprehensive copulas.
The algorithm presented in this section is naturally conceived for one-parameter copulas, but it can be extended to p-parameter copulas, p ≥ 2 (just think of Student's t copula); in this case, p − 1 higher-order correlations or co-moments need to be assigned along with the usual linear correlation in order to calibrate all the parameters.

Extension to multivariate context
The extension of the proposed procedure to dimension d > 2 -finding a joint distributions with d assigned margins and d(d − 1)∕2 distinct pairwise correlations-is not straightforward at all, but this is due to theoretical rather than computational reasons. To explain it, we will refer to a counterexample described in Bergsma and Rudas (2002) and Chaganty and Joe (2006). Let X, Y and Z be three correlated binary rvs, each with support {0, 1} , such that P(X = 1) = P(Y = 1) = P(Z = 1) = 0.5 (nothing changes if we select {1, 2} as common support). Based on the bivariate correlation bounds (see point 3 in Sect. 2) the three correlation coefficients XY , XZ , YZ can lie in the interval [−1, +1] . However, if we choose XY = YZ = 0.4 and XZ = −0.4 , then a trivariate distribution for (X, Y, Z) does not exist. So, a first type of problem is related to the feasibility of the (assigned) correlation matrix P for the discrete d-variate random vector: even if all the bivariate correlation bounds are respected and even if the matrix P reckoning all the pairwise correlations is a valid correlation matrix (a symmetric matrix with all ones on the main diagonal, which is positive semidefinite), nevertheless it may be impossible to construct a random vector with the d assigned margins and correlation matrix P. A second type of problem is related to the nearly lack of copulas able to calibrate-through their parameters-the values of the resulting pairwise correlation coefficients. The typical generalizations of the Frank and Plackett copulas, which we discussed in Sect. 3.2, are still characterized by a unique scalar parameter. Assigning arbitrary-though feasible-values to the pairwise correlations of the discrete random vector, generally would lead to no solution in terms of the dependence parameter. A way to overcome this issue is resorting to a copula whose number of parameters is at least equal to the number of distinct pairwise correlations: an obvious candidate is the Gauss copula; a richer option is represented by the t copula. Alternatively, one can resort to pair copula construction through vines, which are graphical models that represent high-dimensional distributions and can model a rich variety of dependencies. Another limitation of the algorithm of Sect. 3.3 is that it only handles rvs with finite supports: extensions to count variables need to seek for an accurate approximation of the correlation coefficient at step 5, possibly entailing a truncation of the support of the joint distribution of steps 3 and 4, in order to compute the correlation coefficient through a double finite summation.

Pseudo-random simulation
Simulating samples from a bivariate rv with assigned point-scale margins and correlation, built according the procedure described in Sects. 3.2 and 3.3, is straightforward. One can resort to the following general algorithm for meta-copula distributions: 1. Simulate a random sample (u 1 , u 2 ) from the copula C(u 1 , u 2 ; ) , where is the value of the copula parameter recovered through the algorithm of Section 3.3; 2. Set x 1 = F −1 1 (u 1 ) and x 2 = F −1 2 (u 2 ) , where F −1 1 and F −1 2 are the generalized inverse functions of F 1 and F 2 , respectively; 3. (x 1 , x 2 ) is a random sample from the target bivariate distribution, with copula C(u 1 , u 2 ; ) and margins F 1 and F 2 .
Inducing a desired value of correlation between two point-scale… Alternatively, since both X 1 and X 2 have a finite support space, one can resort to the "inversion algorithm" described in (Devroye (1986), pp.85,559) and used in Lee (1997): in this case, one directly considers the joint probability mass function p(i, j) of Eq.
(3) and proceeds as follows: 1. Set N = m 1 × m 2 ; let t , t = 1, … , N be the joint probabilities p ij arranged in descending order; 2. Let y t be the corresponding possible values of Y = (X 1 , X 2 ) , arranged similarly; 3. Define z 0 , z 1 , … , z N in the following way: 4. Simulate a random number u from a standard uniform rv; 5. Return y t , where z t−1 < u ≤ z t .

Application to CUB random variables
A CUB rv X is defined as the mixture, with weights and 1 − , of a shifted binomial with parameters m and and a discrete uniform distribution over the support Corduas (2011) proposed using the Plackett copula in order to construct a one-parameter bivariate distribution from CUB margins; this proposal was later investigated by Andreis and Ferrari (2013), also in a multivariate direction. Here, we reprise and extend these attempts of constructing a bivariate CUB rv, by resorting to the results discussed in Section 3.2. Let us suppose that we want to build a bivariate model with margins X 1 ∼ CUB(m 1 = 5, 1 = 0.4, 1 = 0.8) and X 2 ∼ CUB(m 2 = 5, 2 = 0.7, 2 = 0.3) ; we can find the values of the attainable correlations by using the function corrcheck in GenOrd (Barbiero and Ferrari 2015a), which returns the values min = −0.952003 and max = 0.8640543 . We proceed and select a desired feasible value of correlation between the two CUB variates, say = 0.6. Afterward, we recover the values of Ga (for the Gauss copula), (for the Frank copula), and (for the Plackett copula), according to the iterative procedure illustrated in the previous section. By setting = 10 −7 , we obtain Ga = 0.6898959 , = 5.453455 , and = 11.30106 . Table 3 reports the detailed iterations of the algorithm for the three copula-based models. Note that since > 0 , for the Frank copula we selected a value (1) larger than zero (tentatively, (1) = 1 ) and for the Inducing a desired value of correlation between two point-scale… Plackett copula a value (1) larger than one (tentatively, (1) = 2 ). For the Gauss copula, we set (1) = = 0.6. The three joint pmfs, sharing the same value of linear correlation, are reported in Table 4. It is easy to notice the differences among them. For example, the probability p 23 = P(X 1 = 2, X 2 = 3) takes the values 0.0922, 0.0948, and 0.1008, across the three joint distributions. Figure 2 displays the relationship between the copula parameter (on the x axis) and the corresponding Pearson's correlation of the example bivariate CUB model with Gauss, Frank, and Plackett copula, constructed by applying steps 3 ÷ 5 of the proposed algorithm over a dense grid of uniformly spaced values of . The almost linear trend for the Gauss copula can be easily noted (in this case, the copula parameter is itself a correlation coefficient!), whereas for the other two it is obviously (highly) nonlinear, due also to the unbounded domain for and . From the figure, one can state that X 1 ,X 2 is a monotone concave function of the copula parameter for the Plackett copula, and for the Frank copula only when > 0 : this explains the increasing nature of the sequences (t) in Tables 3b and 3c.

Inferential aspects
If we have a bivariate ordinal sample (x 1t , x 2t ) , t = 1, … , n , that we assume has been drawn from a joint cdf F(x 1 , x 2 ; , 1 , 2 ) = C(F 1 (x 1 ; 1 ), F 2 (x 2 ; 2 ); ) , then parameters' estimation can be carried out through different inferential techniques. First, let us define the log-likelihood function as ( 1 , 2 , ) = ( 1 , 2 , |(x 11 , x 21 ), … , (x 1n , x 2n )) = Inducing a desired value of correlation between two point-scale… with p being the joint pmf of Eq.(3) and n ij the absolute joint frequency of (x i , x j ) , i = 1, … , m 1 , j = 1, … , m 2 . We now present three possible ways to perform point estimation of the distribution's parameters: the first one is the customary maximum likelihood method; the second one is a modification thereof, whose use with copula models is however quite consolidated; the third one is directly suggested by the simulation procedure of Sect. 3.3 and can be regarded to as a by-product of the methodology presented in Sect. 3.

Full maximum likelihood
The most standard estimation technique consists of maximizing with respect to all the three parameters (or parameter vectors) simultaneously. This task can be usually done just numerically (i.e. no closed-form expressions are available for the parameter estimates), by resorting to customary optimization routines. The resulting maximum likelihood (ML) estimates can be thus derived as where is the parameter space.

Two-step maximum likelihood
This technique aims at reducing the computational burden of the previous one, by splitting the original maximization problem into two subsequent (sets of) maximizations in lower dimensions. In the first step, one estimates 1 and 2 separately, by resorting to maximum likelihood estimation, as if the two univariate components of the bivariate sample were independent, i.e. maximizing separately their marginal log-likelihood functions: with n i⋅ and n ⋅j being the observed marginal frequencies of X 1 and X 2 , respectively, thus finding the estimates ̂T S 1 and ̂T S 2 (the superscript "TS" standing for "two-step"). Then, one sets 1 =̂T S 1 and 2 =̂T S 2 in (4) and maximize it with respect to , finding ̂T S . This technique was introduced in a more general context and exhaustively described in Joe and Xu (1996), where it is also named "inference function for margins" (IFM). The authors compared the efficiency of the IFM with the ML by simulation and found that the ratio of the mean square errors of the IFM estimator to the full ML is close to 1. Theoretically, the ML estimator should be asymptotically the most efficient, since it attains the minimum asymptotic variance bound. However, for finite samples, Patton (2006) found that the IFM was often even more efficient than the ML. As a result, IFM is the main estimation method employed in estimating copula models.

Two-step maximum likelihood + method of moment
This method is directly suggested by the algorithmic procedure of Sect. 3.3. First, one estimates the marginal parameters 1 and 2 from sample data x 1i and x 2i , i = 1, … , n , independently maximizing 1 and 2 with respect to 1 and 2 , as for the previous technique, obtaining ̂T S 1 and ̂T S 2 . Then, one considers the maximum likelihood estimates of the marginal cdfs, F j (⋅) = F j (⋅;̂T S j ) , j = 1, 2 , and obtain the estimate of the dependence parameter via the method of moment, by inverting the relationship between the and Pearson's correlation: ̂T SM = g −1 (̂X 1 X 2 ;F 1 ,F 2 ) , where ̂X 1 ,X 2 is Pearson's sample correlation coefficient and TSM stands for "twostep-moment method". The evaluation of g −1 at ̂X 1 X 2 , given F 1 and F 2 , is carried out through the algorithm of Sect. 3.3.
We remark that these three estimation methods are just possible alternatives to be employed for the specific context of copula-based discrete distributions considered in this work. When dealing with copula-based distributions in the continuous case, a straightforward estimation method for the parameter of a specific bivariate parametric copula is the method of moment. It consists in considering a rank correlation between the two rvs (say, Kendall's or Spearman's ), looking for a theoretical relationship between it and the parameter of the copula, and substituting the empirical value of the rank correlation into this relationship to derive an estimate of the copula parameter. The main advantage of this method is that it does not require any assumption about the marginal distributions, since rank correlations are margin-free, i.e. they depend on the copula only and are not affected by the marginal distributions.
When dealing with factor copulas, another consolidated estimation method is represented by a sort of simulated method of moments (SMM), where the "moments" that are used in estimation are not the usual ones, but functions of rank statistics. The SMM estimator is derived as the minimizer of the distance between data dependence measures and dependence measures obtained through Monte Carlo simulation of the model (Oh and Patton 2017). In the context of hierarchical Archimedean copulas, Okhrin and Tetereva (2017) investigate a clustering estimator based on Kendall's by means of Monte Carlo simulations; it is shown to be competitive in terms of statistical properties (bias and variance) and to be computationally advantageous.
However, using these methods would not be convenient in our context, since we know rank correlations lose their nice properties, holding in the continuous set-up, due to the presence of tied values (see, for example, Nešlehová (2007)).

Monte Carlo study
The relative performance of the estimators derived through the three methods described in the previous section, expressed in terms of some statistical indicators such as bias or mean-squared error, can be assessed for finite sample size only via Monte Carlo (MC) simulation. Usually, the estimators of the marginal parameters 1 3 Inducing a desired value of correlation between two point-scale… methods have a very close statistical behaviour; on the contrary, differences are expected to arise among the estimators of the dependence parameter. Here we will examine the joint behaviour and performance of all the parameters' estimators.
For the multivariate case, we recall that the bias of an estimator ̂= (̂1, … ,̂p) � of a p-dimensional vector parameter = ( 1 , … , p ) � is defined as: is said to be an unbiased estimator of if (̂) = for any ∈ . A multivariate generalization of the mean-squared error (MSE) is provided by the MSE matrix (see, for example, Mittelhammer 2013, p.377): The MSE matrix can be decomposed into variance and bias components, analogous to the scalar case. Specifically, MSE is equal to the sum of the covariance matrix of ̂ and the outer product of the bias vector of ̂: The trace of the MSE matrix defines the expected squared distance (ESD) of the vector estimator ̂ from the vector estimand and is equal to: where " tr " denotes the trace of a matrix. Being a scalar, the ESD allows direct and easy comparison among different estimators of the same parameter vector: the lower the value of ESD, the better the estimator. Generally, for the estimator vectors corresponding to the three methods of Sect. 5, such quantities cannot be derived analytically, but an approximation can be obtained through the corresponding MC means computed over S simulation runs: with and ̂t being the value of the vector estimator for run t, t = 1, … , S . The larger the value S, the more accurate the approximation. bias(̂) ∶= (̂) − ; A MC study is designed to assess the relative statistical performance of the three types of vector estimators for a bivariate CUB model under an array of artificial scenarios, which are realized by varying the dependence structure, the CUB marginal parameters, and the sample size. We point out that this study will not allow general conclusions, but is intended merely to demonstrate how the different inferential methods work and check the potential of the proposal. As possible dependence structures, we evaluate the Gauss, Frank, and Plackett copulas: the first with parameter Ga equal to −0.6 or +0.6 ; the second with parameter = −5 or = +5 ; the last with parameter = 1∕4 or = 4 . As CUB marginal parameters, we consider m 1 = m 2 = 5 and m 1 = m 2 = 7 , combined with the following values for the marginal parameter vector M = ( 1 , 1 , 2 , 2 ) � : (0.4, 0.8, 0.4, 0.8) � , (0.4, 0.8, 0.7, 0.3) � , and (0.7, 0.3, 0.7, 0.3) � . For all the 2 × 2 × 3 = 12 combinations above, the sample sizes 50 and 100 are investigated, for a total of 24 artificial settings for each type of copula.
Note that the two-step procedures (Sects. 5.2 and 5.3) require that both empirical marginal distributions assume at least 4 distinct values in order the MLEs to be computed (Piccolo 2003). To ensure feasibility of these two techniques for any simulation run, when they occurred, we simply discarded the samples that do not respect such condition from the MC simulation and kept fixed at S = 1000 the number of feasible samples to draw and to use for the statistical analysis under each artificial setting.
Simulation results reporting the values of Ê SD MC are displayed in Tables 5a, b, and c. Other summary indices computed over the 1, 000 simulation runs, such as the MC mean of the bias vector, are available on request.
For the models with the Gauss copula, the full MLE performs the best in terms of ESD; for the other two estimators, the values of ESD are very close to each other for any setting, pointing out a slight preference towards the two-step MLE for n = 100 . There is also a relevant difference in the values of ESD across the settings; as expected, ESD decreases moving from n = 50 to n = 100 , holding the other parameters fixed. The setting with = (0.7, 0.3, 0.7, 0.3, −0.6) � and n = 100 minimizes the value of the index for all the three estimators; the setting with = (0.4, 0.8, 0.4, 0.8, −0.6) � maximizes it. For the models with the Frank copula, the two-step MLE surprisingly overtakes the full MLE for any setting; the two-step moment estimator shows the worst performance, even if for n = 100 this difference is attenuated. Significant improvement of all the three estimators occur when moving from n = 50 to n = 100 ; the values of distribution's parameters also affect the values of ESD: however, the effect of the model parameters cannot be easily extracted from the results on the examined scenarios.
For the models with the Plackett copula, we have the following interesting result: the MLE is the best performer (smallest Ê SD MC ) when = 1∕4 (corresponding to negative dependence); for the complementary scenarios, corresponding to = 4 and then positive correlation, the two-step MLE is the best performer. The twostep moment estimator has a far worse behaviour than its competitors, apart from one scenario, where it is the second best after full MLE; this is especially apparent for = 4 . Note that in this case, differently from the previous models, the values 1 3 Inducing a desired value of correlation between two point-scale… of Ê SD MC sensibly change moving from = 4 to = 1∕4 , holding fixed the other parameters; this is due to the different magnitude of the two values of here considered, whose choice was related to the different meaning of the copula parameter's values: for the Frank and Gauss copulas, changing the sign of means changing the sign to correlation while keeping its intensity fixed; this does not occur for the Plackett copula, whose parameter has to range within ℝ + (see Fig.2), leading to negative dependence if failing in (0, 1) and to positive dependence if larger than 1. Although in this MC study, the results in terms of statistical performance of the inferential technique proposed in Sect. 5.3 are overall worse than those of the two more consolidated techniques described in Sects. 5.1 and 5.2, nevertheless the former can be still useful for providing starting values for the copula parameter to the maximization routines of the two latter two.
We also measured, on a selection of artificial settings, the total computational times in minutes required by each of the three estimation methods over the S = 1000 MC runs, which are displayed in Table 6. Although the magnitude of the computation time depends on the selected dependence structure, it can be noted that in each setting the two-step method of moment is by far the least time-consuming, followed ML full maximum likelihood method, TS two-step maximum likelihood method; TSM= two-step maximum likelihood method+method of moment Table 6 Overall computation times in minutes for the three estimation methods for some artificial scenarios considered in the MC study Here, n = 100 ; m 1 = m 2 = 7 ; 1 = 0.4, 1 = 0.8, 2 = 0.7, 2 = 0.3 ; " + " and "−" indicates the copula parameter value inducing positive and negative dependence, respectively (see Tables 5a, 5b,

3
Inducing a desired value of correlation between two point-scale… by the two-step maximum likelihood method, and the full maximum likelihood method. We remark that the suggested estimation procedure is faster when moving from the Gaussian copula to the other two dependence structures.

Empirical analysis
In this section, we consider an application of the inferential techniques of Sect. 5 to real data, specifically, the survey data coming from the 2000 International Social Survey Programme (ISSP), which addressed the topic of attitudes to environmental protection and preferred government measures for environmental protection. The prefmod package (Hatzinger and Dittirch 2012) comprises the raw data structured as a dataset with 1, 595 complete observations (one for each respondent) on 11 variables, namely five socio-demographical variables (gender, location of residence, age, country, and education) and six items (with a 5-point rating scale, i.e Likert type).

Respondents from Austria and Great
Britain were asked about their perception of environmental dangers; the questions concerned air pollution caused by cars (variable CAR), air pollution caused by industry (IND), pesticides and chemicals used in farming (FARM), pollution of country's rivers, lakes, and streams (WATER), a rise in the world's temperature (TEMP), and modifying the genes of certain crops (GENE). The answers were given on a 5-point rating scale, with response categories: (1) extremely dangerous, (2) very dangerous, (3) somewhat dangerous, (4) not very dangerous, and (5) not dangerous at all for the environment. We focus on respondents from Austria only and on WATER and GENE items. The joint and marginal empirical distributions of the two items are reported in Table 7. By considering the ratings as numerical values, we can treat this distribution as a sample from a bivariate discrete rv; its sample correlation is equal to 0.2988. The nature of the data suggests using CUB as parametric family for the marginal distributions; the Gaussian, Frank and Plackett copulas are assumed as dependence structures for the joint distribution. Table 8 summarizes the estimation results obtained by applying the full maximum likelihood method and both the two-step methods. From Table 8, it is easy to note that within each dependence structure, the three estimation methods provide estimate values for the marginal parameters which are very close to each other; the estimates of the dependence parameter are slightly different. For example, under the Gaussian dependence structure, the values of the estimates of the dependence parameter Ga range between 0.343 and 0.361by the way, this confirms a positive and moderate dependence between the two observed variables, as suggested by the value of the sample correlation. Differences in terms of maximized log-likelihood emerge across the three dependence structures here examined. Among the three models, the bivariate CUB with Plackett copula (and parameters set equal to the corresponding MLEs) is that providing the best fit ( = −2006.317).
Alternative ways of comparing the goodness of fit of the three models can be taken up; for example, one can consider Aitchison distance between the empirical and theoretical probabilities = ∑ ∑ (log(p ij ∕p ij ) −L) 2 , where L = ∑ ∑ log(p ij ∕p ij )∕k , with k being the total number of points of the support of the bivariate rv. The minimum value of the Aitchison distance is achieved by the model with the Plackett copula ( = 3.830 ); the model with the Frank copula provides = 4.089 ; the Gauss copula returns = 4.390 . Other possible distances or divergences between discrete distributions are mentioned in Fossaluza et al. (2018).
An evaluation in absolute terms of the goodness-of-fit of this bivariate model can be carried out by resorting to the chi-squared statistic defined as 2 = ∑ ∑ (n obs ij − n theo ij ) 2 ∕n theo ij , with n obs ij and n theo ij indicating the observed and the theoretical frequencies, respectively. In order to accomplish the requirement that each n theo ij has to be not smaller than 5, we can proceed and collapse the last two ordered categories (4 and 5) for each variable, thus obtaining the new observed joint distribution in Table 9. There, we also reported the theoretical joint frequencies corresponding to the "best" bivariate CUB model, for which the Chi-squared statistic above takes the value 21.68, with a p-value 0.0168 (the test statistic, under the null hypothesis that the data come from the bivariate CUB model, asymptotically follows a chi-squared distribution with 16 − 1 − 5 = 10 degrees of freedom). This means that the goodness of fit of the model is hardly satisfactory; at the significance level 1%, we do not reject the null hypothesis. Looking at Table 9, discrepancies between

3
Inducing a desired value of correlation between two point-scale… observed and theoretical joint (and also marginal) frequencies are visible to the unaided eye. One can also consider not using a parametric model for the two margins and estimating them nonparametrically. In the opposite direction, other families of bivariate distributions may be tested for fit improvement; for example, cumulative link models could be used for the univariate margins (Agresti and Kateri 2019). Reasonably, introducing covariates (gender, education, age, location of residence) for the marginal (and dependence) parameters would likely increase the fit. However, we remark that the aim of this section was not so much to evaluate the fitting of bivariate (CUB) models to real data as to illustrate and compare different estimation techniques, one of which derived through a correlation-matching procedure. Besides, we remark that pooling cells, although being a viable alternative to obtain accurate p-values in some instances, should be preferably performed before the analysis is made in order to obtain a statistic with the appropriate asymptotic reference distribution; otherwise, it may distort the purpose of the analysis (Maydeu-Olivares and García-Forero 2010). Alternatively, one could resort to resampling methods (e.g. bootstrap), but unfortunately, existing evidence suggests that resampling methods do not yield accurate p-values for the 2 statistic (Tollenaar and Mooijart 2003).

Conclusions
In this work, we showed how to build a joint probability distribution with assigned discrete point-scale margins enjoying a target (feasible) value of correlation. We proposed a two-step copula-based approach: first, one selects a copula function and constructs a bivariate distribution preserving the assigned margins; then, one adjust the value of the copula parameter in order to achieve the target correlation: this leads to the elaboration of an iterative procedure whose accuracy can be set a priori, differently from other approaches in the literature that are based on some rearrangement of very huge samples drawn independently from the two margins.
The approach is designed to work with any one-parameter copula family, provided that it encompasses the entire range of dependence, thus allowing the use of other dependence structures than the Gaussian, whose use in the social sciences has been overriding. Although in this paper we considered three exchangeable and radially symmetric copulas, this feature is not necessary for the algorithm functioning. Moreover, being comprehensive is a property that makes the algorithm applicable to any feasible correlation value; however, if one is only concerned, for example, with positive correlations, then he/she can use non-comprehensive copulas as well, such as Gumbel or Clayton.
As said, the algorithm is specifically conceived for one-parameter copulas, but it can be extended to two or more-parameter copulas: in this case, some higher comoments need to be assigned along with linear correlation in order to calibrate the additional parameters.
The extension of the proposed procedure to dimension d > 2 , discussed in Sect. 3.4, is not straightforward, being limited by theoretical rather than computational reasons, related to the properties of Pearson's correlation matrix for a d-variate random vector.
Future research will explore these aspects.

Supplementary material
The R code implementing the algorithm of Section 3, along with the example of Section 4, is provided as supplementary material here: https:// tinyu rl. com/ ASTB-D-19-00215.
Funding Open access funding provided by Università degli Studi di Milano within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.