* This example do-file tests if the "price" variable from 
* the "auto" dataset distributed with Stata is power law distributed
* the methods implemented come from:
* A. Clauset, C.R. Shalizi, and M.E.J. Newman, 
* "Power-law distributions in empirical data" SIAM Review 51(4), 661-703 (2009)
* (arXiv:0706.1062, doi:10.1137/070710111)

* The methods implemented here assume that the variable under study
* is continuous; the methods for discrete power law model are not implemented

* Notice that this do-file follows closely the above paper's methos
* and it will be hardly possible to understand the do-file 
* without the close reading of the paper

* I have tested the programs used in this do file, but
* they still may contain errors; the help files are not ready yet. 

sysuse auto,clear
*cd C:\Users\Michal\Desktop\Desktop\wealthLists\0programs\programs_clean_versions

* setting a lower limit for searching for x0
qui sum price
local xmin=r(min)

* estimates the lower bound x0 on power-law behaviour
* see details in Clauset et al. (2009)
plx0 price, x0start(`xmin')
local ntail=r(ntail)
local x0=r(x0)

* fits a power law model with a given x0
plfit price, x0(`x0')
local alpha = e(alpha_cont)
local pl_lik = e(likelihood)

* calculated log-likelihoods for power law model 
* needed for Vuong's model selection tests performed below
qui gen double pl_ll = -`alpha'*ln(price) + log(`alpha'-1) + (`alpha'-1)*ln(`x0') ///
if price>=`x0'

qui cou
local obs = r(N)
* plots data with a power law model fitted
plplot price, obs(`obs') alpha(`alpha') tailobs(`ntail')

local nreps = 199
* performs the parametric bootstrap goodness of fit test for power law model
plgof price, x0(`x0') nreps(`nreps')

* fits log-normal model by ML to be compared below with power law model
* if the model does not coverge, other optimization techniques may be used
* see: help ml (maximize_options, expecially "difficult" and "technique" options)
lognfit2 price if price>=`x0',x0(`x0') difficult
mat re = e(b)
local v = re[1,2]
local m = re[1,1]
* generate log-likelihods for log-normal model to be used in Vuong's test below
qui gen double logn_ll = -ln(price) - ln(sqrt(2*_pi)) - ln(`v') 	///
		- .5*(`v'^(-2))*( ln(price) - `m')^2 - ln(1-normal((ln(`x0') -`m')/`v')) ///
		if price>=`x0'

* fits exponential model by ML to be compared below with power law model
exponfit price if price>=`x0',x0(`x0')
mat res=e(b)
local lam=res[1,1]
* generate log-likelihoods for exponential model 
qui gen double exp_ll = ln(`lam') - (`lam' * price) + `lam'*`x0' ///
		if price>=`x0'

* fits stretched exponential model by ML to be compared with power law below
stretchexponfit price if price>=`x0',x0(`x0')
mat res=e(b)
local lambda=res[1,1]
local beta = res[1,2]
* generate log-likelihoods for stretched exponential model
qui gen double strexp_ll=ln(`lambda') - `lambda' * price^`beta' + ////
	`lambda'*(`x0')^`beta' + ln(`beta') + (`beta' - 1)*ln(price) ///
	if price>=`x0'

* fits power law with exponential cut-off model	
plcutfit price if price>=`x0',x0(`x0') from(`alpha' `lam', copy)
mat res=e(b)
local alpha=res[1,1]
local lambda = res[1,2]
local plcut_lik=e(ll)
* computes incomplete gamma function needed to compute log likelihoods for power law with exp. cut-off
incompletegamma `alpha' `lambda'
local gamma = r(gamma)
* generate log-likelihoods for the power law with exp. cut-off
qui gen double plcut_ll = ///
	-`alpha'*ln(price) - `lambda'*price + (1-`alpha')*ln(`lambda') ///
	-ln(`gamma') if price>=`x0'

* Vuong's model selection tests

* 1) power law vs log-normal
vuong, ll1(pl_ll) ll2(logn_ll)
* 2) power law vs exponential
vuong, ll1(pl_ll) ll2(exp_ll)
* 3) power law vs stretched exponential
vuong, ll1(pl_ll) ll2(strexp_ll)

* 4) power law vs power law with exponential cut-off
* assumes that densities are correctly specified

loc lr = - 2 * (`pl_lik' - `plcut_lik')
di chi2tail(1, `lr')