We propose a novel approach for the computation of the probability distribution of a counting variable linked to a particular kind of hierarchical multivariate copula function called a clusterized homogeneous copula. Here, the problem considered is very complex in a high-dimensional setting. As is common practice for large-dimensional problems, we restrict ourselves to positive orthant dependence and we define that copula on clusterized data, allowing us to reduce the dimension of the problem. This approach approximates a multivariate distribution function of heterogenous variables with a distribution of a fixed number of homogeneous clusters, organized through a clustering method as proposed in a 2011 paper by the authors. To compute the probability density function of the counting variable linked to a clusterized homogeneous copula, we propose an algorithm, implemented in Matlab code. We compare this probability density function with that computed through the Panjer recursion approach and the limiting Gaussian and Archimedean approaches, which are commonly used for high-dimensional problems. The scalability of the algorithm is also studied. As an application, we study the problem of evaluating the distribution of losses related to the default of various types of counterparty in a credit risk exposed portfolio.