EconBase
← Back to paper

Composite distributions in the social sciences: A comparative empirical study of firms' sales distribution for France, Germany, Italy, Japan, South Korea, and Spain

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

77,158 characters · 13 sections · 61 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Composite distributions in the social sciences: A comparative empirical study of firms' sales distribution for France, Germany, Italy, Japan, South Korea, and Spain

\titlerunning{Composite distributions for firm sizes}

\authorrunning{Ramos et al.}

abstractWe study 17 different statistical distributions for sizes obtained from the classical and recent literature to describe a relevant variable in the social sciences and Economics, namely the firms' sales distribution in six countries over an ample period. We find that the best results are obtained with mixtures of lognormal (LN), loglogistic (LL), and log Student's $t$ (LSt) distributions. The single lognormal, in turn, is strongly not selected. We then find that the whole firm size distribution is better described by a mixture, and there exist subgroups of firms. Depending on the method of measurement, the best fitting distribution cannot be defined by a single one, but as a mixture of at least three distributions or even four or five. We assess a full sample analysis, an in-sample and out-of-sample analysis, and a doubly truncated sample analysis. We also provide the formulation of the preferred models as solutions of the Fokker--Planck or forward Kolmogorov equation. \ \\ JEL Codes: C61, D39, L25.

Introduction

The study of the size variables in the social sciences and Economics has a long tradition like the use of a power-law in the study of the income distribution Par1896, the application to the upper tail of the size distribution of cities Sin36,zipf1949human, see e.g., also the recent article by Par21. On its side, the distribution of firms has been the subject of early studies as the classic one by R. Gibrat Gib31, where the lognormal distribution has been proposed, based on the now called Gibrat's law or the law of proportionate effect. Classical references have proposed the lognormal distribution for firm size and variations SimBon58,Qua66,Cla79,CabMat03. In the city size distribution, the lognormal has been proposed, e.g., by Par85 and eeckhout2004gsl.

More recent publications dealing with firm sizes' distribution and related matters could be, for example, Ree01,Ree02,Ree03,reed2004double, GioLevRan11,GuoXuCheWan13,IoaSko13,Tod17,CorMorPer17,KwoNad19,BanChiPrePueRam19, Su19,GuaTos18,GuaTos19,GuaTos19b,PueRamSan20,PueRamSanArr20, CamRam21,Reg21,Ish21. Recent articles on the issue about whether the upper tail is lognormal or Pareto could be, e.g., Per05,ClaShaNew09,BeeRicSch11,BeeRicSch13,Bee15,BeeRicSch17,ChuDicNad19,SchTre19,Bee22,ChaFla22. Also, there are studies aiming at the explanation of economically important quantities by mixtures of distributions BelCleBul14,Kun20, being quantities for which sub-populations of it may exist McLPee00. Let us mention as well that this research has roots on previous one of several years ago GonRamSan14,Pue15,RamSan15,Ram17.

For the theoretical explanations of the so-called Zipf's law and Gibrat's law standard references are Gab99,Gab09 and Gib31, and in modern terms, they can be understood, respectively, as distributions that are solutions to Fokker--Planck or forward Kolmogorov equations associated to It\^{o} differential equations Ord74,Gar04,ItoMcK96,Kyp06 with, respectively, reflective lower barrier Har85 or not, but both with constant diffusion and drift terms. Let us remark that the solutions to Fokker--Planck or forward Kolmogorov equations are mathematical concepts, but have been accepted by the Economics community after Gab99,Gab09 as plausible economical explanations of the observed empirical laws, as those papers are widely cited afterwards.

However, there is a more recent reference on the subject that sheds light on a fact known after Dupire Dup93,Dup94 but somehow overlooked in the literature that many probability density functions can be obtained in that way. In fact, following PenPueRamSan22, let $y=\ln(x)$ be the logarithm of the firm's sales variable $x>0$. We assume that its evolution or dynamics is governed by the It\^{o} differential equation (see, e.g., Ord74,Gar04)

equation*[equation* omitted — 60 chars of source]

where $B_{t}$ is a standard Brownian motion (Wiener process) (see, e.g., ItoMcK96,Kyp06 and references therein). The quantity $a(y_{t},t)$ corresponds to the diffusion term, and $b(y_{t},t)$ to the drift term. This process can be associated to the forward Kolmogorov equation or Fokker-Planck equation for the time-dependent probability density function (PDF) (conditional on the initial data) $f(y,t)$ (see also Gab99,Gab09):

equation[equation omitted — 170 chars of source]

Since the probability density function $f(y,t)$ is evolving on time and perhaps there is no limiting stationary distribution, let us propose a way of solving the Fokker--Planck equation for the mentioned $f(y,t)$ by specifying the diffusion term and the drift term (see, e.g., Otu21 for another recent approach to time-dependent solutions of the Fokker--Planck equation). In fact, if we take $a(g,t)=s^2$, where $s>0$ is a real constant, then by choosing

equation[equation omitted — 154 chars of source]

where ${\rm cdf}(y,t)$ is the cumulative distribution function (CDF) corresponding to $f(y,t)$, this solves the Fokker--Planck equation for $f(y,t)$ Dup93,Dup94. However, we remark that we do not claim that such a solution with this choice of $a(y,t)$ and $b(y,t)$ is unique, only that $f(y,t)$ is a solution by construction. Also, $b(y,t)$ might have bounded discontinuities in the variable $y$, in a finite number of points in the domain GikSko07. And a third remark is that we may add to the expression of $b(y,t)$ above a term of the form $h(t)/f(y,t)$, where $h(t)$ is an arbitrary function of $t$.

Thus we intend to fill the gap from theory (that many PDFs can be generated in this way) to empirical facts (what probability density function or functions is or are observed empirically), employing this empirical work, for the manufacturing firm's sales distribution for six countries in an ample period.

We offer a panoramic view of the subject for manufacturing firm's sales distribution since we have at hand ample databases for several countries, covering the full sales distribution, and not only the upper tail. For that, we select some distributions that have been proposed recently for size studies, including firms' sizes, that satisfy that the PDFs are differentiable at least three times in the PDF parameters, and that have a CDF expressible using elementary or special functions known to date.

The rest of the paper is organized as follows. Section (ref) describes the distributions to be compared. Section (ref) describes the databases used. Section (ref) offers the results. Finally, we end with some conclusions.

The distributions

As mentioned earlier, $x>0$ will denote the sales of the firm in question. The first distribution that we will consider, as a baseline model, is the usual lognormal distribution (LN), given by

equation[equation omitted — 143 chars of source]

where $\mu\in\mathbb{R}$, $\sigma>0$. The corresponding CDF is

equation*[equation* omitted — 123 chars of source]

where ${\rm erf}$ denotes the error function associated to the standard normal distribution.

The second distribution in our study is the double Pareto lognormal (DPLN) Ree02,Ree03,reed2004double,Man09,giesen2010size, given by

eqnarray[eqnarray omitted — 390 chars of source]

where $\mu\in\mathbb{R}$, $\alpha,\beta,\sigma>0$ are the four parameters of the distribution. It has the property that it approximates different power-laws in each of its two tails: $f_{\rm DPLN}(x)\approx x^{-\alpha}$ when $x\to\infty$ and $f_{\rm DPLN}(x)\approx x^{\beta}$ when $x\to 0$. The body is approximately lognormal, although it is not possible to delineate exactly the switch between the lognormal and the power-law behaviors since the DPLN distribution is the log version of the convolution of an asymmetric double Laplace with a normal distribution.

The third distribution in our study is the Generalized Beta of the second kind (GB2) McD84,McDXu95,KleKot03 and depends on four parameters with probability density function:

equation[equation omitted — 130 chars of source]

where $B(p,q)=\int_0^1 t^{p-1}(1-t)^{q-1}\,dt$ is the Beta function, and $a,b,p,q>0$ are the four distribution parameters.

A way of modelling firm size distribution with semi-nonparametric densities is found in CorMorPer17. The log-semi-nonparametric (LNSNP) density is defined using the Hermite polynomials until degree 4 and the lognormal density as follows:

eqnarray[eqnarray omitted — 171 chars of source]

where $\mu\in\mathbb{R}$, $\sigma>0$, $d_1,d_2,d_3,d_4\in\mathbb{R}$ are the distribution parameters, $z$ is a shorthand for ${\displaystyle \frac{\ln(x)-\mu}{\sigma}}$ and

eqnarray[eqnarray omitted — 123 chars of source]

are the cited Hermite polynomials.

The $\ell$-mixtures of lognormal distributions ($\ell$LN) for $\ell\geq 2$ are defined, for example, as follows McLPee00,GuaTos18,GuaTos19b,GuaTos19,KwoNad19,Su19, PueRamSan20,PueRamSanArr20,CamRam21

eqnarray[eqnarray omitted — 259 chars of source]

where $0\leq p_1,\dots,p_{\ell-1}\leq 1$, $0\leq 1-\sum_{j=1}^{\ell-1}p_j\leq 1$, and $\mu_i\in\mathbb{R}$, $\sigma_i>0$, $i=1,\dots,\ell$.

The specific lognormal mixtures that we will consider are

eqnarray[eqnarray omitted — 173 chars of source]

Likewise, the loglogistic (LL) distribution has the probability density function $$ f_{{\rm LL}}(x;\mu,\sigma) =\frac{\exp{\left(-\frac{\ln(x)-\mu}{\sigma}\right)}} {x\sigma \left(1+\exp{\left(-\frac{\ln(x)-\mu} {\sigma}\right)}\right)^2} $$ where $\mu\in\mathbb{R},\sigma>0$ and $x>0$. The $\ell$-mixtures of loglogistic distributions ($\ell$LL) for $\ell\geq 2$ can be defined, analogously, as McLPee00,PueRamSan20

eqnarray[eqnarray omitted — 259 chars of source]

where $0\leq p_1,\dots,p_{\ell-1}\leq 1$, $0\leq 1-\sum_{j=1}^{\ell-1}p_j\leq 1$, and $\mu_i\in\mathbb{R}$, $\sigma_i>0$, $i=1,\dots,\ell$.

The specific loglogistic mixtures that we will consider are

eqnarray[eqnarray omitted — 173 chars of source]

Some other mixtures could be considered. We adapt three other mixtures that appeared firstly in MasPueRam20 and afterwards in MasRam21, called 2St12, 2St39, and 3St, to yield the 2LSt12, 2LSt39, and 3LSt, and other supplementary two that will be called 4LSt, 5LSt. They are as follows. If the log version (to correctly deal with a size variable as it is firms' sales) of the non-stan\-dar\-di\-zed Student's $t$ distribution is JohKotBal95

eqnarray[eqnarray omitted — 245 chars of source]

where $\mu\in\mathbb{R},\sigma>0$ and $\nu>0$ is the number of degrees of freedom parameter, and $\Gamma(\cdot)$ denotes the Gamma function. This distribution, without logs, has been used to study the log-growth rates of the size of German cities SchTre16.

The $\ell$-mixtures of LSt distributions for $\ell\geq 2$ can be defined, analogously, as McLPee00

eqnarray[eqnarray omitted — 312 chars of source]

where $0\leq p_1,\dots,p_{\ell-1}\leq 1$, $0\leq 1-\sum_{j=1}^{\ell-1}p_j\leq 1$, and $\mu_i\in\mathbb{R}$, $\sigma_i>0$, $\nu_i>0$, $i=1,\dots,\ell$.

The specific log-Student's $t$ mixtures that we will consider are

eqnarray[eqnarray omitted — 466 chars of source]

Note that we fix the number of degrees of freedom a priori. There are several reasons for that. First, the maximum likelihood estimation that we will perform afterwards becomes much more stable than with free degrees of freedom parameters for $\ell\geq 2$. The choice a priori of the parameters of the degrees of freedom in the mixtures allows us to avoid convergence problems of the numerical estimation process. Second, with these choices, we break down the existing symmetry in the estimation parameters (concerning a scenario with free parameters of degrees of freedom in which the interchanging of the parameters of the components lead to the same distribution) so identification is achieved in this case. And third, the versions of these distributions without logs have worked very well when studying log-growth rates of city sizes MasPueRam20, and the log-returns of stock indices worldwide MasRam21.

The data sets

We use data of sales of manufacturing firms (measured in constant thousands of USA dollars, the reference year being 2005) of France (FR), Germany (DE), Italy (IT), Japan (JP), South Korea (KR) and Spain (ES), from the ORBIS database. For France, we have data for the years 2005-2014, for Germany, we have 2010-2014, for Italy we have 2006-2011 and for Japan, South Korea, and Spain we have 2005-2014. ORBIS is the largest database of firms' financial data in the world. We used the 2016 edition, which still contains financial data for more than 100 million firms worldwide. It is characterized not only by the size and completeness of the data but also by the simultaneous analysis of multiple countries because the data of each country is described in the same format. While the latest data from the 2016 edition of ORBIS is 2015, the number is limited because they are still being collected. ORBIS contains data for approximately 10 years. In this paper, we focus on France (FR), Germany (DE), Italy (IT), Japan (JP), South Korea (KR), and Spain (ES) as countries with sufficient data for this period.

We will use the data in three different ways in our empirical application. First, the whole data without cut-offs put us in the worst situation to obtain a good fit. Second, we will split the data into two (approximately) equally distributed subsamples: the 75% of each of the data sets will be the in-sample data, and the other 25% of each sample will be the out-of-sample data. We will assess the fit of the out-of-sample data into the estimated in-sample distributions. Third, for each full sample we will drop the 10% of the bottom observations (lower tail) and the 0.1% of the top observations (upper tail) to improve the measures of skewness and kurtosis of the log data as it is suggested in Fio20, and more importantly, to show that the tails of the full and doubly truncated data are very difficult to model directly. We also consider all the manufacturing firms without separating by type of industry, following again the idea to put us in the worst situation to obtain a good fit, since then the sample size will be maximum and the goodness-of-fit using standard statistics will be increasingly demanding.

The descriptive statistics for the full samples can be seen in Table (ref). It can be observed that the log data presents skewness and kurtosis different from that of the normal distribution, so it is guessed that a single lognormal will not be a good fit for the sales variable for these data sets.

table[table omitted — 5,477 chars of source]

We show likewise the descriptive statistics of the in-sample and out-of-sample data sets in Tables (ref) and (ref), respectively. We can see similar descriptive statistics between these two subsamples and with regards to the full samples as well.

table[table omitted — 5,474 chars of source]
table[table omitted — 5,454 chars of source]

And finally, we show the descriptive statistics of the doubly truncated (tt) data sets in the way explained before, in Table (ref). Now the descriptive statistics do vary somehow, in particular, the skewness and kurtosis of the log data are slightly closer to the ones of the normal distribution, as suggested in Fio20.

table[table omitted — 5,318 chars of source]

Empirical results

We will summarize in this Section the results of fitting different models to each different type of data, separating them into different Subsections. But there are aspects of the method followed that are common to all the analyses. First, we have only taken into account distributions that could have the regularity conditions of, for example, Kie78,BasMcL85,NewMcF94,CasBer02,AtiGarMunVil07 and avoiding distributions whose probability density is not three times continuously differentiable in the parameters and distributions that have not a CDF expressible in terms of elementary or special functions known to date. Also, the maximum likelihood estimation (MLE) that we have carried out has been done using the command {\tt mle} of {\sc MATLAB}$^\circledR$, which relies on the Nelder--Mead simplex algorithm NelMea65. We have chosen this to treat all the distributions on an equal footing, instead of choosing an EM algorithm for the mixture models McLPee00 and using the above-cited command {\tt mle} for the other distributions. Independently of the estimations obtained by the {\tt mle} command, we have checked that the estimations correspond to a maximum of the log-likelihood by varying the initial starting points. We have checked as well that the maxima are so around intervals centred on the estimate and with a width of eight standard errors. The standard errors are computed following the procedure of EfrHin78 and McCVin03.

To assess the goodness-of-fit of the different models to describe the data, we have computed standard Kolmogorov--Smirnov (KS), Cr\'amer--von Mises (CM), and Anderson--Darling (AD) statistics. Let us recall how these statistics are computed.

itemize• The Kolmogorov-Smirnov (KS) statistic KolmogorovAN1933 compares the empirical distribution function ${\rm cdf}_n(x)=\frac{1}{n}\sum_{i=1}^n\mathds{1}_{(0,x]}(x_i)$ with the distribution function of the fitted model ${\rm cdf}_{B}(x;\hat{\theta})$ and is given by the maximal deviance \begin{equation*} KS=\sup_{x\in (0,\infty)}|{\rm cdf}_n(x)-{\rm cdf}_B(x;\hat{\theta})| \end{equation*} and $n$ is the sample size. • The Cram\'er--von Mises (CM) statistic Cra28,Mis28 is defined by \begin{align*} CM&=n\int_{0}^{\infty}({\rm cdf}_n(x)-{\rm cdf}_{B}(x;\hat{\theta}))^2\mathrm{d}{\rm cdf}_{B}(x;\hat{\theta})\\ &=\frac{1}{12n} +\sum_{i=1}^n\left(\frac{2i-1}{2n}-{\rm cdf}_{B}(x_{(i)};\hat{\theta})\right)^2, \end{align*} • The Anderson-Darling (AD) statistic AndDar54 is defined by \begin{align*} AD&=n\int_{0}^{\infty}\frac{({\rm cdf}_n(x)-{\rm cdf}_{B}(x;\hat{\theta}))^2}{{\rm cdf}_{B}(x;\hat{\theta})(1-{\rm cdf}_{B}(x;\hat{\theta}))}\mathrm{d}{\rm cdf}_{B}(x;\hat{\theta})\\ &=-n-\sum_{i=1}^n\frac{2i-1}{n}\left(\log {\rm cdf}_{B}(x_{(i)};\hat{\theta})+\log\left(1-{\rm cdf}_{B}(x_{(n+1-i)};\hat{\theta})\right)\right), \end{align*} where $x_{(1)}\leq\ldots\leq x_{(n)}$ are the observed ordered data.

The smaller the value of these statistics or distances, the most favoured the evaluated model.

In order to select models, we resort to standard information criteria very well adapted to the ML estimation. They are:

itemize• The Akaike Information Criterion (AIC) Aka74,BurAnd02,BurAnd04, defined as $$ AIC=2k-2\ln L^* $$ where $k$ is the number of parameters of the distribution and $\ln L^*$ is the corresponding (maximum) log-likelihood. The minimum value of AIC corresponds (asymptotically) to the minimum value of the Kullback--Leibler divergence, so a model with the lowest AIC is selected from among the competitors. • The Bayesian or Schwarz Information Criterion (BIC) Sch78,BurAnd02,BurAnd04, defined as $$ BIC=k \ln(n)-2 \ln L^* $$ where $k$ is the number of parameters of the distribution, $n$ the sample size, and $\ln L^*$ is as before. The BIC penalizes more heavily the number of parameters used than does the AIC. The model with the lowest BIC is selected according to this criterion. • The Hannan--Quinn Information Criterion (HQC) HanQui79,BurAnd02,BurAnd04, defined as $$HQC=2 k \ln(\ln(n))-2 \ln L^* $$ where $k$ is the number of parameters of the distribution, $n$ the sample size, and $\ln L^*$ is as before. The HQC implements an intermediate penalization of the number of parameters when compared to the AIC and BIC. The model with the lowest HQC is selected according to this criterion.

Full samples' results

We start showing the empirical results for the full samples. The used distributions can be estimated almost always, with the following exceptions: The GB2 and LNSNP for JP 2005, the 5LN, 3LL, and 4LL for IT 2007-2011, the 4LL for JP 2014, the 5LL for IT 2007-2010 and JP 2014, the 3LSt for JP 2010, and the 5LSt for IT 2007-2011. For the rest of the cases, all estimations correspond to a well-behaved maximum with corresponding standard errors that are shown on different sheets of a supplementary Excel book.\footnote{This file is available at \href{https://doi.org/10.7910/DVN/PHHJM3} {https://doi.org/10.7910/DVN/PHHJM3}.}

Next, the computation of the goodness-of-fit statistics KS, CM, and AD yield the summarized outcomes in Table (ref).

table[table omitted — 1,704 chars of source]

From Table (ref) we can see that several distributions never correspond to the minimum of either KS, CM, or AD, amongst them, the lognormal distribution, and several others. We will consider the 51 possibilities times three statistics each, which amounts to 153 cases in total, and also in what follows. The top-six of minimum statistics models are the 5LL (52 times out of 153), 5LSt (27), 5LN (23), 4LL (15), 4LSt (13), 3LN (11), showing that the mixtures yield the best results and with an increasing number of components in the mixture. This is reasonable, as with more components the fit is supposed to be better. This is our next task to see how these models are selected according to information criteria, that penalize the number of parameters.

We show in Table (ref) the summarized results for the full samples' case. We can see that the five-top selected models are the 4LSt (45 times out of 153), the 5LN (31), the 3LN (23), the 4LN (18), and the 5LSt (13). The non-mixture models are never selected, but some mixtures like the 2LL, 3LL, 2LSt12, 2LSt39, not either. One could argue that there is a slight tendency to have better and better information criteria with a higher number of components. We conjecture that this is due to the difficulty of modelling the tails of the distribution of the full data sets, as it happened for example in Ram22 and references therein when studying the upper tail of world billionaires' data. We will consider the doubly truncated data sets in the way explained before in Subsection (ref) and then we will see the differences in the selected models concerning those in this Subsection.

table[table omitted — 1,710 chars of source]

A study of in-sample and out-of-sample assessment

In this subsection, we study the descriptive power of different models when we separate each full sample into two subsamples: in-sample (75% of the full data) and out-of-sample (25% of the full data), the separation being made in an identically distributed way. We have checked whether the two sub-samples so obtained come from the same distribution with standard KS, CM, and AD tests, and the global result is that at the 5% level the null hypothesis of coming from the same distribution is not rejected 143 out of 153 times, so it becomes a very reasonable assumption. The details are shown on a sheet of the cited supplementary Excel book.\footnote{Again, this file is available at \href{https://doi.org/10.7910/DVN/PHHJM3} {https://doi.org/10.7910/DVN/PHHJM3}.}

Then, we have fitted the 17 distributions used in this paper to the in-sample data sets and computed the standard errors as before. The cases that we have not been able to estimate are the following: The GB2 and LNSNP for JP 2005, the 4LN for IT 2007-2010, the 5LN for IT 2007-2011, the 3LL for IT 2007-2011, the 4LL for IT 2006 and 2009 and JP 2014, the 5LL for IT 2007-2011 and JP 2010, 2013-2014, and the 5LSt for IT 2007-2011.

We have computed the KS, CM, and AD statistics of the out-of-sample data sets when evaluated on the in-sample estimated distributions, and also the corresponding log-likelihoods, AIC, BIC, and HQC information criteria. It is known that this procedure shows more parsimony in the selected models than the full sample estimations.

table[table omitted — 1,720 chars of source]

We can see in Table (ref) that there is not a clear ranking in the classification of the different distributions with regards to the outcome of KS, CM, AD statistics when considering non-mixture or mixture distributions, and amongst these last ones, there is not a clear ranking between a lower number of components and higher number of components, as the lowest statistics are scattered across different distributions.

As with regards to the information criteria, we show in Table (ref) the corresponding summary. We see that the top-five minimum information criteria are obtained for the 3LN (43 times out of 153), the 2LN (21), the 2LL (20), the 4LSt (16), and the 4LN (13). Thus we observe that none of the non-mixture models is amongst the top minimum information criteria distributions and that amongst the mixture models, there is a majority of minimum information criteria with only three (3LN) or two (2LN, 2LL) components, shown in this way parsimony in the number of parameters of the descriptive power of out-of-samples of different distributions.

table[table omitted — 1,715 chars of source]

Doubly truncated data sets and distributions

In this subsection, we will consider the doubly truncated data sets and distributions. That is, we drop from the full data sets the lowest 10% of observations and the top 0.1% of the observations, as is suggested for example in Fio20. Also, the distributions are modified as follows JohKotBal94. Suppose we have a random variable $X$ that is distributed according to some PDF $f(x)$, with CDF ${\mathrm{cdf}}(x)$, both of which have support equal to $(0,\infty)$. Suppose we wish to know the PDF of the random variable after restricting the support to $[a,b]$ with $0<a<b<\infty$. Then

equation[equation omitted — 100 chars of source]

where $g(x)=f(x)I(\{a\leq x\leq b\})$ and $I$ is the indicator function. If the support of the non-truncated distribution is $(-\infty,\infty)$, the procedure is entirely similar. The formulae for the CDFs of our considered distributions are available upon request. For the non-mixture models we perform this truncation as just explained. For the mixture models we perform the truncation of each of the components first, and after that taking the mixture of them. There is the alternative of taking first the mixture of the components and after that, doubly truncating the so obtained mixture. We feel that the first alternative provides better results as the estimations are more stable and more of them out of the possible ones are obtained. In all cases, the doubly truncated distributions corresponding to the ones described in Section (ref) will be appended with \lq\lq tt\rq\rq.

table[table omitted — 1,062 chars of source]

We observe in Table (ref) that the top-four lowest statistics correspond to the 5LLtt (63 times out of 153), 5LSttt (26), 4LLtt (25), and 5LNtt (11) and that the non-mixture distributions never provide the lowest statistics. Thus there is a tendency to have lower statistics with five-component mixture distributions (5LLtt and 5LSttt) followed by a four-component one (4LLtt), although next is also a five-component one (5LNtt). This is not unexpected for the reasons we have pointed out before, namely that the goodness-of-fit does not penalize the higher number of parameters.

To consider that, let us have a look at the information criteria. We show in Table (ref) the summary of the results. We see that the top-five most selected distributions are the 5LSttt (28 times out of 153), 3LNtt (14), 2LLtt (12), 5LLtt (11), and 3LSttt (8). The non-mixture distributions are selected 5 times in total out of 153, but never the LNtt. Thus we do not obtain a clear parsimony in this case with the number of the parameters, as there is an ample selection of five-component mixture distributions. This may be caused by the difficulty of modelling the upper tail of the distributions, even after having removed the top 0.1% of the observations, in line with Ram22 and references therein.

table[table omitted — 1,083 chars of source]

Fokker--Planck equations for the preferred models

We will briefly show the solution for $b(y,t)$ of ((ref)) associated to the Fokker--Planck equation ((ref)) in the case of the 4LSt (full samples) and 5LSttt (doubly truncated samples). The corresponding preferred model for the in-sample out-of-sample analysis, namely the 3LN, can be treated in a similar way and it is already presented in PenPueRamSan22 (with a variable $g$ instead of $y$, it is easy to perform the substitution).

Let us begin with the 4LSt. We denote first the natural logarithm of $x>0$ by $y=\ln(x)\in(-\infty,\infty)$. The PDF of the ordinary Student's $t$ distribution in the variable $y$ is given by $$ f_{{\rm St}}(y;\mu,\sigma,\nu)= \frac{\Gamma\left(\frac{\nu+1}{2} \right)}{\Gamma\left(\frac{\nu}{2}\right)\sqrt{\pi \nu}\sigma} \left(1+\frac{1}{\nu}\left(\frac{y-\mu}{\sigma} \right)^2 \right)^{-\frac{\nu+1}{2}}\nonumber $$ and the corresponding CDF is given by $$ {\rm cdf}_{{\rm St}}(y;\mu,\sigma,\nu)= \frac{1}{2}+ \frac{\Gamma\left(\frac{1+\nu}{2}\right)} {\Gamma\left(\frac{\nu}{2}\right)\sqrt{\pi\nu}\sigma} (y-\mu)\,{}_2F_{1}\left(\frac{1}{2},\frac{1+\nu}{2}, \frac{3}{2},-\frac{(y-\mu)^2}{\nu \sigma^2} \right) $$ where ${}_2F_{1}$ denotes the Gaussian or ordinary hypergeometric function. Let us denote for the sake of brevity the corresponding time-dependent density function by

eqnarray[eqnarray omitted — 291 chars of source]

Also, let us denote for convenience the following expression: $$ k(y;\mu,\sigma,\nu,s)=\dot{\mu} +(y-\mu)\frac{\dot{\sigma}}{\sigma} -\frac{s^2(1+\nu)(y-\mu)}{2((y-\mu)^2+\nu \sigma^2)} $$ where $\mu\in\mathbb{R}$ and $\sigma,\nu>0$. $\mu$ and $\sigma$ are supposed to depend smoothly on $t$ (explicit dependence is omitted for notational simplicity) and $\nu$ is fixed (a priori), the dot means derivative with respect to $t$, and $s>0$ is a real constant. In this expression, taking $\nu\rightarrow\infty$ yields the corresponding expression for the normal distribution, as it should be, presented in PenPueRamSan22 (for the variable $g$ instead of $y$). Let us define also the quantities

eqnarray[eqnarray omitted — 417 chars of source]

and the time-dependent posterior probabilities

eqnarray[eqnarray omitted — 393 chars of source]

Then, if we select $a(y,t)=s^2$, with $s>0$ being a real constant, the $b(y,t)$ given by ((ref)) turns into

eqnarray[eqnarray omitted — 284 chars of source]

so that we obtain that $f(y,t)=j_{\rm 4St}(y,t)$ is in this case a solution of the corresponding Fokker--Planck equation ((ref)). The sign of the drift term is indefinite on this occasion. This is related to the results of CamRam21,PenPueRamSan22.

For the doubly truncated 5LSt (5LSttt) case, let us denote again first $y=\ln(x)$. Then, we doubly truncate the corresponding Student's $t$ distribution to give ($y_{\rm min}=\ln(x_{\rm min})$ and $y_{\rm max}=\ln(x_{\rm max})$) $$ f(y;y_{\rm min},y_{\rm max},\mu,\sigma,\nu)= \frac{f_{{\rm St}}(y;\mu,\sigma,\nu) I(\{y_{\rm min}\leq y\leq y_{\rm max}\})}{{\rm cdf}_{{\rm St}}(y_{\rm max};\mu,\sigma,\nu)-{\rm cdf}_{{\rm St}}(y_{\rm min};\mu,\sigma,\nu)} $$ with support $-\infty<y_{\rm min}\leq y\leq y_{\rm max}<\infty$ and the corresponding CDF being given by $$ {\rm cdf}(y;y_{\rm min},y_{\rm max},\mu,\sigma,\nu)= \frac{({\rm cdf}_{{\rm St}}(y;\mu,\sigma,\nu)-{\rm cdf}_{{\rm St}}(y_{\rm min};\mu,\sigma,\nu))I(\{y_{\rm min}\leq y\leq y_{\rm max}\})} {{\rm cdf}_{{\rm St}}(y_{\rm max};\mu,\sigma,\nu)-{\rm cdf}_{{\rm St}}(y_{\rm min};\mu,\sigma,\nu)} $$ The corresponding time-dependent mixture is then

eqnarray[eqnarray omitted — 458 chars of source]

Then, let us denote again for convenience the following expression: $$ k(y;y_{\rm min},y_{\rm max},\mu,\sigma,\nu,s)= \frac{s^2}{2f}\frac{\partial f}{\partial y} -\frac{1}{f}\left(\dot{\sigma}\frac{\partial{\rm cdf}}{\partial \sigma}+\dot{\mu}\frac{\partial{\rm cdf}}{\partial \mu}+\dot{y}_{\rm max}\frac{\partial{\rm cdf}}{\partial y_{\rm max}}+\dot{y}_{\rm min}\frac{\partial{\rm cdf}}{\partial y_{\rm min}}\right) $$ where $\mu,y_{\rm min},y_{\rm max}\in\mathbb{R}$, $\sigma>0$, and are supposed to depend smoothly on $t$ (explicit dependence is omitted for notational simplicity), and $\nu>0$ is kept fixed; the dot means derivative with respect to $t$, and $s>0$ is a real constant. Let us define also the quantities

eqnarray[eqnarray omitted — 673 chars of source]

and the time-dependent posterior probabilities

eqnarray[eqnarray omitted — 537 chars of source]

Then, if we select $a(y,t)=s^2$, being again $s>0$ a constant, the $b(y,t)$ given by ((ref)) turns into

eqnarray[eqnarray omitted — 508 chars of source]

so that we obtain that $f(y,t)=j_{\rm{5Sttt}}(y,t)$ is in this case a solution of the corresponding Fokker--Planck equation ((ref)). The sign of the drift term is indefinite also on this occasion. This development for mixtures of doubly truncated distributions is new as far as we know. Also, it is easy to substitute the Student's $t$ in these mixtures by another distribution to give the corresponding expressions for other mixtures defined on $(-\infty,\infty)$ or doubly truncated ones.

These solutions of the Fokker--Plank equation for the preferred time-dependent density functions might lead to a new way of simulating and/or perhaps forecasting the firm's sales of the different countries, leaving that research to a subsequent paper.

Conclusions

We have considered ORBIS data sets for six different countries, on ample time intervals, and a total of 51 samples, so we believe that the results are quite robust.

We have studied the fit of 17 different statistical models obtained from the classic and recent literature on sizes distributions in the social sciences and Economics, and applied them to the firm size distribution measured by sales, with mixtures of distributions or not. Out of them, there is a classical model, the lognormal (LN), that is strongly not selected always for all samples. When the full data sets are considered, the most often selected distribution is the 4LSt with 11 parameters, and it is not a distribution with the highest number of parameters out of that considered (for example, the 5LN, 5LL, and 5LSt have all 14 parameters). When performing in-sample out-of-sample analysis, the descriptive power of the 3LN has the best results out of our 17 distributions, and again the LN is never selected. Finally, when we doubly truncate the data and the distributions, the 5LSttt is the most often selected distribution, and never the LNtt. This may be due to the that there may be firms which operate by different mechanisms and then conform to some sub-populations in each of the samples, as in other contexts BelCleBul14,Kun20 and as described in McLPee00. Depending on the method of measurement, for the best-fit distribution, it is necessary to consider a mixture of at least three distributions or even four or five ones.

Thus we have two clear results: First, in none of the situations the lognormal or variations thereof is the selected distribution, and depending on the situation one mixture distribution KwoNad19,BanChiPrePueRam19,GuaTos19,GuaTos19b,Su19,PueRamSan20,PueRamSanArr20 is preferred to other ones, when describing firm size distribution measured by sales. Second, it is clearly stated the difficulty of modelling the tails of the full or doubly truncated data sets: A higher number of components in the mixture should be taken into account to obtain better information criteria for these samples.

Also, the mixtures considered in this paper, when taken as time-dependent ones, can be shown to be solutions of a Fokker--Planck equation with constant diffusion term as described in the introduction of this paper, following the lines of CamRam21,PenPueRamSan22. We defer the use of such newly introduced solutions in simulating or forecasting to another paper.

Author contributions

Shouji Fujimoto: Conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, software, supervision, validation, visualization, writing-original draft. Atushi Ishikawa: Conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, software, supervision, validation, visualization, writing-original draft. Till Massing: Conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, software, supervision, validation, visualization, writing-original draft. Takayuki Mizuno: Conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, software, supervision, validation, visualization, writing-original draft. Arturo Ramos: Conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, resources, software, validation, visualization, writing-original draft.

Competing interests statement

The authors declare to have no competing interests concerning the research carried out in this article.

Data availability statement

The Edition of 2016 ORBIS database used in this study has been obtained upon payment from Bureau van Dijk, and we have signed a confidentiality agreement so we cannot disclose the data sets. For the sake of replication, interested researchers can buy the same edition of such database.

Acknowledgments

We thank Profs. Josefina and Rafael Cabeza-Laguna for their help in the computations for the estimation of the doubly truncated mixture 5LSttt.