/* This program helps test whether a given data set fits a
power law distribution. To begin, you must already have estimate a plfit
for your data and computed a KS test for your data (e.g. ksmirnov mag = freq);
you will input both the alpha from your estimation and your observations into
the simulation. After the simulation, you will need the KS statistic for your
estimate of alpha.
This program comes with NO WARRANTY.
(C) Miles Townes
mdtownes@gwmail.gwu.edu
version 0.1 - 3 Dec 2008
*/
*program drop pltest //comment out this line for first run
program define pltest, rclass
version 10.1
syntax [, obs(integer 1) xmin(real 1) alpha(real 1)]
drop _all
tempvar mag c freq
set obs `obs'
gen double `c' = (`alpha'-1)*(`xmin'^(`alpha'-1))
gen double `mag' = `xmin'*((1-runiform())^(-1/(`alpha'-1)))
gen double `freq' = (`c'/(`alpha'-1))*(`mag'^-(`alpha'-1))
//you can use this same formula on real data to generate
//a cdf distribution
ksmirnov `mag' = `freq'
gen double ks = r(D_1)
end
simulate ks, reps(2500) saving(kstest3) : pltest, obs(381) xmin(1000) alpha(1.49913)
//reps() should be 2500 for p-value accurate to .01
//saving() identifies the name of the new dataset
//obs() should be equal to your dataset observations
//xmin() determines the lower bound of the dataset
//alpha is the slope parameter given by your estimate of the power-law fit
count if _sim_1>.9952
//replace .9952 with your KS statistic as appropriate.
//Your p-value is the result of the count divided by
//the number of reps (e.g. 2500). If the number is below
//your threshold - e.g. .05 - you can reject the power-law
//distribution.