A simple test for power-law behavior

Abstract

In this article, I propose a new test for power-law behavior. The statistical test, pwlaw, is locally optimal if the possible alternative distributions are contained in the Pareto type (IV) family. After deriving the test, I examine four classical datasets: the frequency of unique words in an English text (Moby Dick); the human populations of U.S. cities; the frequency of U.S. family names; and the peak gamma-ray intensity of solar flares. I show that in the first case there is no indication of any power-law behavior and that in the second and fourth cases there is evidence in that regard.

1 Introduction

A continuous random variable is said to have a power-law (PWL) behavior if it has a probability density of the form

$g(x)=\frac{\alpha }{\mu }{\left(\frac{x}{\mu }\right)}^{-(\alpha +1)}$

1

where α > 0 is the shape parameter that determines the tail of the distribution and µ > 0 is the location parameter where the “heavy tail” starts. ¹ The law given above is due to Pareto (1897), but in the particular case of α = 1, it is occasionally named after Zipf (1949).

The surveys by Brock (1999), Mitzenmacher (2004), and Newman (2005) make clear that the assumption of PWL behavior has been routinely made over the years by many authors in a number of disciplines. Also, in most of that literature, the use of PWL distributions is typically defended by simply pointing out the approximate linear nature of a doubly logarithmic plot or by reporting a simple regression. However, as noted by, for example, Clauset, Shalizi, and Newman (2009) and Urzúa (2011), those arguments are flawed.

Fortunately, there are also in the literature better procedures to test for PWL behavior. The most common, which can be traced back at least to Quandt (1964), is to estimate, for a given dataset, not only the Paretian law but also some competing alternative distributions. Once those estimations are made, the resulting fits are ranked according to some statistical criteria. For instance, in their rather influential article on the empirics of PWL behavior, Clauset, Shalizi, and Newman (2009) suggest as alternatives to that distribution four others: the lognormal, the exponential, the stretched exponential, and the PWL with cutoff.

Other testing procedures have been proposed as well. In particular, using a methodology that is quite common in the econometrics literature, Urzúa (2000) presents a test for Zipf’s law that is derived after considering a more general distribution that nests the density given in (1). Afterward, a locally optimal test for Zipf’s law is derived by means of Rao’s (1948) score test, also known in econometrics as the Lagrange multiplier test. ²

The more general density posited in that article is of the form

$f(x)=\frac{\alpha }{\mu }{\left(1+\frac{x-\mu }{\sigma }\right)}^{-(\alpha +1)}x>\mu$

2

where, aside from the two parameters introduced earlier, σ > 0 is the scale parameter. This distribution is known in the literature by several names, but here it will be simply called Pareto (II), following the terminology of Arnold (see the reference in Arnold [2015]); likewise, Pareto’s basic law, given in (1), will be called from now on Pareto (I).

Using the Pareto (II) distribution, Urzúa (2000) derives a test for Zipf’s law by testing the statistical hypothesis H ₀ : σ = µ, α = 1. Rao’s score test is used because it requires estimation only under the null hypothesis. Note that if the null is reduced to H ₀ : σ = µ, then, as Goerlich (2013) shows, one can even produce a test for the Pareto (I) distribution. Nevertheless, that test would not be appropriate in the case of testing for PWL behavior, because, as I will show in the next section, the densities (1) and (2.1) have a similar heavy-tail behavior.

2 A test for PWL distributions

Following again the terminology in Arnold (2015), this article assumes that the alternatives to the Pareto (I) distribution belong more generally to the Pareto (IV) family, whose density is given by

$f(x)=\frac{\alpha }{{\gamma }^{\sigma }}{\left\{1+{\left(\frac{x-\mu }{\sigma }\right)}^{1/\gamma }\right\}}^{-(\alpha +1)}{\left(\frac{x-\mu }{\sigma }\right)}^{1/\gamma -1}x>\mu$

3

where γ > 0 is the inequality parameter. Note that (2.1) is obtained from (2.1) when σ = µ and γ = 1 and that the Pareto (IV) distribution function is given by

$F(x)=1-{\left\{1+{\left(\frac{x-\mu }{\sigma }\right)}^{1/\gamma }\right\}}^{-\alpha }x>\mu$

For purposes of this article, the parameter µ, which is present in the three densities given above, is assumed to be fixed by statistical design. This assumption is quite reasonable from a practical point of view, aside from the fact that it is theoretically required. ³ Thus, after defining w = x − µ, we can rewrite the Pareto (IV) density in a more compact form as

$F(w)=\frac{\alpha }{{\lambda }^{\sigma }}{\left\{1+{\left(\frac{w}{\sigma }\right)}^{1/\gamma }\right\}}^{-(\alpha +1)}{\left(\frac{w}{\alpha }\right)}^{1/\gamma -1}w>0$

4

Note that the tail behavior of a Pareto (IV) distribution is determined not only by α, as is the case for the Pareto (I) and (II) distributions, but also by the new parameter γ. This is so because, using L’Hôpital’s rule and (2),

$\begin{array}{c}\underset{w\to \infty }{\lim }\frac{1-F(w)}{w^{-\alpha /\gamma }}=\underset{w\to \infty }{\lim }\frac{\gamma f(w)}{\alpha w^{-\alpha /\gamma -1}}\\ =\underset{w\to \infty }{\lim }{\left\{1+{\left(\frac{w}{\sigma }\right)}^{1/\gamma }\right\}}^{-(\sigma +1)}{\left\{{\left(\frac{\sigma }{w}\right)}^{1/\gamma }\right\}}^{-(\sigma +1)}{\sigma }^{\alpha /\gamma }\\ ={\sigma }^{\alpha /\gamma }\end{array}$

5

Thus, the Pareto (IV) family can be approximated, on the right tail, by w^−α/γ . This implies that, for a given PWL parameter α, the tail of the distribution can become as short as desired by reducing the value of γ from one to some small but positive value. Likewise, for a given α, the tail can be made even heavier by setting the value of γ greater than one.

Let θ = (θ ₁ , θ ₂ , θ ₃) = (σ, γ, α) be the vector of parameters, and let ${\left\{w_i\right\}}_i^n=\text{1}$ be a random sample. Then, the corresponding log-likelihood function is

$\mathcal{l}(\theta )=n\log \left(\frac{\alpha }{\gamma \sigma }\right)-(\alpha +1){\displaystyle \sum _{i=1}^n\log }\left\{1+{\left(\frac{w_i}{\sigma }\right)}^{1/\gamma }\right\}+(1/\gamma -1){\displaystyle \sum _{i=1}^n\log }\left(\frac{w_i}{\sigma }\right)$

For later use, it is necessary to find its gradient and the expected value of its Hessian. The elements of the gradient are given by

$\frac{\partial \mathcal{l}}{\partial \sigma }=-\frac{n}{\gamma \sigma }+\frac{\alpha +1}{\gamma \sigma }{\displaystyle \sum _{i=1}^n\frac{{(w_i/\sigma )}^{1/\gamma }}{1+{(w_i/\sigma )}^{1/\gamma }}}$

6

$\frac{\partial \mathcal{l}(\theta )}{\partial \gamma }=-\frac{n}{\gamma }+\frac{1}{{\gamma }^2}{\displaystyle \sum _{i=1}^n\log (w_i/\sigma )+\frac{\alpha +1}{{\gamma }^2}{\displaystyle \sum _{i=1}^n\frac{{(w_i/\sigma )}^{1/\gamma }\log (w_i/\sigma )}{1+{(w_i/\sigma )}^{1/\gamma }}}}$

7

$\frac{\partial \mathcal{l}(\theta )}{\partial \alpha }=\frac{n}{\alpha }-{\displaystyle \sum _{i=1}^n\log \left\{1+{(w_i/\sigma )}^{1/\gamma }\right\}}$

8

The second-order derivatives are also easy to find, but they are omitted here for lack of space. A less obvious task is to calculate their expected values. Fortunately, they have already been derived in the meticulous article by Brazauskas (2002). Thus, the elements of the information matrix

$I_{ij}(\theta )=-E_{\theta }\left\{\frac{{\partial }^2\mathcal{l}(\theta )}{\partial {\theta }_i\hspace{0.17em}\partial {\theta }_j}\right\}i,j=1,2,3$

are also known.

The stage is now set to devise a test for the PWL distribution. This is accomplished here by employing Rao’s score test in the case of the null hypothesis H ₀ : σ = µ, γ = 1. As a first step to that end, we find the restricted maximum-likelihood estimate for the nuisance parameter α by equating (8) to zero and replacing w_i with x_i − µ:

$\widehat{\alpha }={\left\{\frac{1}{n}{\displaystyle \sum _{i=1}^n\log (x_i/\mu )}\right\}}^{-1}$

9

Consequently, the restricted score vector is given by $\widehat{d}=(d_{\text{1}},d_{\text{2}},0{)}^{\prime }$ , where the first two components can be found using (8) in (4)–(5) and imposing σ = µ and γ = 1. After some simple algebra, they are found to be

$d_1=\frac{n\widehat{\alpha }}{\mu }-(\widehat{\alpha }+1){\displaystyle \sum _{i=1}^n\frac{1}{x_i}}$

10

$d_2=-n+\widehat{\alpha }{\displaystyle \sum _{i=1}^n\log \left(\frac{x_i-\mu }{\mu }\right)}-(\widehat{\alpha }+1){\displaystyle \sum _{i=1}^n\frac{\mu }{x_i}\log \left(\frac{x_i-\mu }{\mu }\right)}$

11

Before presenting the corresponding information matrix, I need to introduce some notation. Let $\Gamma (a)={\displaystyle \int { }_0^{\infty }} t^{a-}{ }^{\text{1}}{\text{e}}^{-t}\text{d}t,\hspace{0.17em}\psi (a)={\Gamma }^{\prime }(a)/\Gamma (a)$ , and ψ′(a) = dψ(a)/da be, respectively, the gamma, digamma, and trigamma functions. Also, to have cleaner expressions in what follows, I define p(a) ≡ ψ(a) − ψ(1) − 1 and q(a) ≡ ψ′(a) + ψ′(1). Given that notation, and after imposing σ = µ and γ = 1 in the information matrix for the Pareto (IV) distribution in Brazauskas (2002, 165), we obtain the following restricted information matrix:

$\widehat{I}=n\left\{\begin{array}{ccc}\frac{\widehat{\alpha }}{{\mu }^2(\widehat{\alpha }+2)}& \frac{\widehat{\alpha }p(\widehat{\alpha })+1}{\mu (\widehat{\alpha }+2)}& \frac{1}{\mu (\widehat{\alpha }+1)}\\ \frac{\widehat{\alpha }p(\widehat{\alpha })+1}{\mu (\widehat{\alpha }+2)}& \frac{\widehat{\alpha }\{p{(\widehat{\alpha })}^2+q(\widehat{\alpha })+2(p(\widehat{\alpha })+1)}{\widehat{\alpha }+2}& \frac{p(\widehat{\alpha })}{\widehat{\alpha }+1}\\ -\frac{1}{\mu (\widehat{\alpha }+1)}& \frac{p(\widehat{\alpha })}{\widehat{\alpha }+1}& \frac{1}{{\widehat{\alpha }}^2}\end{array}\right\}$

12

Finally, keeping in mind that $\widehat{d}=(d_{\text{1}},d_{\text{2}},0{)}^{\prime }$ is defined by (10) and (9) and that $\widehat{I}$ is given in (10), we define the score test for PWL distributions as

$\text{PWL}={\widehat{d}}^{\prime }{\widehat{I}}^{-1}\widehat{d}\stackrel{a}{\sim }{\chi }_2^2$

Under the null, PWL is asymptotically distributed as a chi-squared with two degrees of freedom. The statistic can be readily computed using the pwlaw command, with the following syntax:

pwlaw varname [in] [, mu( # )]

Note that the program automatically drops any x_i ≤ µ.

How well does the test statistic behave in the case of small samples? Table 1 presents evidence that the asymptotic critical values can be safely used when 1 ≤ α ≤ 3, the typical range for the shape parameter, and n ≥ 100, as is the case for most applications. Otherwise, the significance points in table 1 can be used as approximations, or for sharper results, the pwlaw command can be used to perform a specific Monte Carlo simulation.

Table 1.

Significance points for PWL

α	Level	n =	20	25	30	40	60	100	∞
1	5%		5.64	5.70	5.77	5.81	5.85	5.93	5.99
	10%		4.19	4.28	4.35	4.40	4.46	4.53	4.61
1.5	5%		5.51	5.60	5.71	5.74	5.84	5.88	5.99
	10%		4.06	4.18	4.25	4.32	4.41	4.48	4.61
2	5%		5.44	5.54	5.63	5.70	5.80	5.86	5.99
	10%		3.96	4.08	4.17	4.26	4.36	4.45	4.61
2.5	5%		5.39	5.49	5.60	5.68	5.79	5.83	5.99
	10%		3.89	4.01	4.10	4.21	4.32	4.41	4.61
3	5%		5.35	5.45	5.55	5.65	5.78	5.85	5.99
	10%		3.83	3.95	4.05	4.16	4.28	4.39	4.61

NOTE: Simulations using 100,000 replications.

3 Applications

We now exemplify the use of the PWL test in the case of four interesting datasets studied by Newman (2005) and Clauset, Shalizi, and Newman (2009). ⁴ They are listed in table 2, together with seven values associated with each of them. Regarding the columns, note that although x _max is clearly the maximum value in each sample, ${\widehat{x}}_{\text{min}}$ is not the minimum value but rather an estimated lower bound on PWL behavior. ⁵ Given the value of that bound, the next smaller number in the sample corresponds to µ (to enforce x > µ), while n _tail is the effective sample size.

Table 2.

Tests of PWL behavior

	n	xmax	${\widehat{x}}_{\text{min}}$	µ	ntail	PWL	p-value
Count of word use	18,855	14,086	7	6	2,958	189	0.00
Population of cities	19,447	8,008,654	52,457	52,360	580	0.16	0.92
Frequency of surnames	2,753	2,502,020	111,919	109,432	239	3.32	0.19
Solar flare intensity	12,773	231,300	323	322	1,711	0.87	0.65

The first dataset in table 2 refers to the number of times that unique words occur in a classical piece of English literature, the novel Moby Dick by Herman Melville. Our statistic can be computed with the command pwlaw x, mu(6). As shown in the table, the value of the statistic thus obtained is quite large in that particular case, and hence the associated p-value is quite small. Thus, according to the PWL test, the null hypothesis of PWL behavior is strongly rejected. This result is at odds with many previous studies on word frequencies, ranging from Zipf (1949) to Clauset, Shalizi, and Newman (2009).

The second case corresponds to the human populations of U.S. cities as recorded by the U.S. Census Bureau in 2000. As shown in the table, and according to the proposed test, there is convincing evidence of PWL behavior. This result is in line with most of the literature on the size of cities. Note also that this finding does not contradict the claim in Urzúa (2000) that U.S. cities do not follow Zipf’s law, because the latter is just an extreme case of a PWL distribution (when α = 1).

The third sample is made of the frequencies of U.S. family names in the 1990 U.S. Census. As reported in the table, there is moderate evidence of PWL behavior, because the null hypothesis would not be rejected at the typical significance levels. A very similar finding is reported by Clauset, Shalizi, and Newman (2009).

Finally, the fourth sample listed in table 2 is made of the peak gamma-ray intensity of solar flares, measured from orbit between 1980 and 1989 (see Newman [2005] for more details). Given the p-value reported in the table, there is good evidence of PWL behavior in this case. Note that a theoretical model that renders precisely that behavior is due to Lu and Hamilton (1991).

4 Concluding remark

The Pareto (IV) distribution used in this article belongs to the Feller–Pareto family (see Arnold [2015]), whose general density is of the form

$f(x)=\frac{\Gamma ({\gamma }_1+{\gamma }_2)}{\gamma \sigma \hspace{0.17em}\Gamma ({\gamma }_1)\Gamma ({\gamma }_2)}{\left\{1+{\left(\frac{x-\mu }{\sigma }\right)}^{1/\gamma }\right\}}^{-({\gamma }_1+{\gamma }_2)}{\left(\frac{x-\mu }{\sigma }\right)}^{({\gamma }_1+{\gamma }_2)-1}x>\mu$

where γ ₂ > 0 and, to use the notation given earlier, γ ₁ = α. Thus, following the same procedure given before, we can derive a new, more general test statistic. But because of its algebraic intricacy, it is not presented here. In any case, the gain from that generalization is not totally apparent, because, as can be checked in the same way as in (3), the Feller–Pareto and Pareto (IV) distributions have similar tail behavior.

Supplemental Material

Supplemental Material, st0610 - A simple test for power-law behavior

Supplemental Material, st0610 for A simple test for power-law behavior by Carlos M. Urzúa in The Stata Journal