/* This program adapts the RANHT method from MATLAB code available at
http://www.santafe.edu/~aaronc/powerlaws/
This is a highly simplified version of the original RANDHT method.
This program comes with NO WARRANTY.
(C) Miles Townes
mdtownes@gwmail.gwu.edu
version 0.1 - 24 Nov 2008
version 2.0 - 3 Dec 2008
*/
program drop plsim //comment out this line for first run
program define plsim, rclass
version 10.1
syntax [, obs(integer 1) xmin(real 1) alpha(real 1)]
drop _all
set obs `obs'
tempvar mag c freq
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))
gen double mag = `mag'
gen double freq = `freq'
end
simulate mag freq, reps(1000) saving(plsim1) : plsim, obs(1) xmin(1000) alpha(1.5)
//reps() determines the number of observations in your dataset
//saving() identifies the name of the new dataset
//obs() should be 1 for all simulations
//xmin() determines the lower bound of the dataset
//alpha() is the slope parameter of the power-law fit
rename _sim_1 mag
rename _sim_2 freq
//simulate returns "x" values as a variable named _sim_1
gen logmag = log10(mag)
gen logfreq = log10(freq)
*gen year = _n //a pseudo-variable for the author's research
*twoway (scatter logmag year) //optional, probably not useful
twoway (scatter logfreq logmag)
//the power-law scatterplot; note that the slope of the graph is alpha - 1