# Power Law Analysis # Exampe Code # Christian Breunig # Oct 15, 2009 ## libraries library(e1071) # for kurtosis library(Lmoments) # for L-kurtosis library(MASS) # for getting the optimal bin with # generate some data x<-c(rweibull(1000, shape=1.1, scale = 2),-rweibull(1000, shape=1.3, scale = 1)) ## histogram # determine bandwidth estimate and optimal number of bins x<-na.omit(x) # drop the missing values h<-width.SJ(x) ob<-(max(x)-min(x))/h # histogram with optimal bins and normal curve #pdf("hist-t.pdf", height = 6, width = 9, paper = "special") hi<-hist(x, col="grey", probability=T, nclass=ob, main="Histogram of Growth Rate in Outlays", xlab = "Growth Rate", ylab = "Probability Density") lines(density(rnorm (100000, mean=mean(x),sd=sd(x)))) #dev.off() ## Plots # category midpoints and densities of the histograms (separate positive and negative tails) himid<-hi$mids hiden<-hi$density nem<-as.double(abs(himid[himid <0 & hiden!=0])) nec<-cumsum(as.double(hiden[himid<0 & hiden!=0])) pom<-as.double(abs(himid[himid>0 & hiden!=0])) poc<-as.double(hiden[himid>0 & hiden!=0]) poc<- sort(cumsum(sort(poc)),decreasing=T) # note the following plots can be prettified by plotting on log scales etc # log-log plot #pdf("llplot.pdf", height = 6, width = 6, paper = "special") plot(log10(c(nem,pom)),log10(c(nec,poc)), type="n",main="",xlab="Log of Growth Rate",ylab="Log Probability",cex.axis=0.9) points(log10(pom),log10(poc), pch=16,col=1) reg<-lm(log10(poc)~log10(pom)) abline(coef(reg),lwd=1,lty=1,col=1) points(log10(nem),log10(nec),pch=1,col=1) reg2<-lm(log10(nec)~log10(nem)) abline(coef(reg2),lwd=1,lty=2,col=1) legend(x="topright",title="Tail",lwd=1,col=1,lty=1:2,legend=c("positive","negative"),cex=0.8,bty="n") reco1<-round(reg$coeff[2],digits=2) rer1<-round(summary(reg)$r.squared,digits=2) mtext(bquote(beta[pos]==.(reco1) ~R^2 == .(rer1)), side=3,line=1) reco2<-round(reg2$coeff[2],digits=2) rer2<-round(summary(reg2)$r.squared,digits=2) mtext(bquote(beta[neg]==.(reco2) ~R^2 == .(rer2)), side=3,line=0) mtext("Log-Log Plot", side=3,line=2.5, cex=1.5) #dev.off() #semi-log plot #pdf("xlplot.pdf", height = 6, width = 6, paper = "special") plot(c(nem,pom),log10(c(nec,poc)), type="n",main="",xlab="Growth Rate",ylab="Log Probability",cex.axis=0.9) points(pom,log10(poc), pch=16,col=1) reg<-lm(log10(poc)~pom) abline(coef(reg),lwd=1,lty=1,col=1) points(nem,log10(nec),pch=1,col=1) reg2<-lm(log10(nec)~nem) abline(coef(reg2),lwd=1,lty=2,col=1) legend(x="topright",title="Tail",lwd=1,col=1,lty=1:2,legend=c("positive","negative"),cex=0.8,bty="n") reco1<-round(reg$coeff[2],digits=2) rer1<-round(summary(reg)$r.squared,digits=2) mtext(bquote(beta[pos]==.(reco1) ~R^2 == .(rer1)), side=3,line=1) reco2<-round(reg2$coeff[2],digits=2) rer2<-round(summary(reg2)$r.squared,digits=2) mtext(bquote(beta[neg]==.(reco2) ~R^2 == .(rer2)), side=3,line=0) mtext("Semi-Log Plot", side=3,line=2.5, cex=1.5) #dev.off()