Skip to content

Instantly share code, notes, and snippets.

@Protonk
Created November 23, 2010 19:09
Show Gist options
  • Save Protonk/712312 to your computer and use it in GitHub Desktop.
Save Protonk/712312 to your computer and use it in GitHub Desktop.
PPP and GDP relationships
#The following code requires a few packages to be installed, namely animation, WDI, and ggplot2.
#You can install them by typing install.packages("ggplot2"), inserting the package name. This
#requires an internet connection, because it is downloading a package from a mirror.
#Loads a library to access World Bank data automatically.
library(WDI);
#Loads a graphics library which makes nicer graphs than base R.
library(ggplot2);
#How does this change across years?
#Names the indicators so I don't need to type them in individually.
wdi.ind<-c("NY.GDP.PCAP.CD","NY.GDP.PCAP.PP.CD","NY.GDP.MKTP.CD");
#Loads up the data from the WDI API over the internet. Requires an internet connection
GDP.giant<-WDI(country="all",indicator=wdi.ind,start=1980,end=2009,extra=TRUE);
#Names the variables to something more tractable
names(GDP.giant)[4:6]<-c("GDP.per","GDP.PPP","GDP");
#Removes Aggregates (region averages).
GDP.giant<-subset(GDP.giant,GDP.giant$region != "Aggregates");
#Builds the dependent variable, lambda.
GDP.giant$lambda<-abs(GDP.giant$GDP.PPP-GDP.giant$GDP.per)/GDP.giant$GDP.per;
#Removes the US and Canada, both of which are pretty uninformative for us.
GDP.giant<-GDP.giant[which(GDP.giant$region != "North America"),];
#Changes region names into categorical variables and renames them slightly.
GDP.giant$region<-as.factor(GDP.giant$region);
levels(GDP.giant$region)<-c("East Asia & Pacific","Europe & Central Asia","Latin America\n& Caribbean","Middle East & N. Africa","South Asia","Sub-Saharan Africa");
#Plots
overall.plot<-qplot(log(GDP.per),lambda,data=GDP.giant,geom="point",alpha=I(0.3))+stat_smooth(method="lm");
byregion.plot<-qplot(log(GDP.per),lambda,data=GDP.giant,colour=region,geom="point",alpha=I(0.3),log="x")+facet_wrap(~region)+scale_y_continuous(limits=c(0,4))+scale_x_continuous('')+opts(aspect.ratio = 0.5);
lambda.year.plot<-qplot(year,lambda,data=GDP.giant,colour=region,geom="smooth")+scale_y_continuous(limits=c(0,2));
#Build a dataframe for the value of Beta Across the years
i<-0
beta.years<-matrix(0,2,30);
#This is a control structure which calculates a linear model for each year in the dataset and fills
#in the point estimate for the slope and the stardard error of that estimate into a matrix which
#we turn into a data frame.
for(i in sort(unique(GDP.giant$year))) {#
beta.years[1,i-min(unique(GDP.giant$year))+1]<-coef(summary(lm(lambda~log(GDP.per),data=subset(GDP.giant,year==i))))[2,1]#
beta.years[2,i-min(unique(GDP.giant$year))+1]<-coef(summary(lm(lambda~log(GDP.per),data=subset(GDP.giant,year==i))))[2,2]
}
GDP.beta<-data.frame(t(beta.years),row.names=c(sort(unique(GDP.giant$year))))
names(GDP.beta)<-c("Beta","S.E.");
#construct confidence intervals for beta from computed standard errors and the (probably wrong) assumption of normally distributed errors.
GDP.beta$lower<-with(GDP.beta,Beta-(qnorm(0.975)*S.E.));
GDP.beta$upper<-with(GDP.beta,Beta+(qnorm(0.975)*S.E.));
#Plot confidence intervals
betaint.plot<-qplot(as.numeric(dimnames(GDP.beta)[[1]]),lower,data=GDP.beta,geom="line",main="Beta and 95% Confidence Intervals")+geom_line(mapping=aes(as.numeric(dimnames(GDP.beta)[[1]]),Beta,data=GDP.beta),size=2,colour="blue")+geom_line(mapping=aes(as.numeric(dimnames(GDP.beta)[[1]]),upper,data=GDP.beta))+scale_y_continuous("Beta")+scale_x_continuous('');
#The below code is adopted from https://github.com/hadley/ggplot2/wiki/Using-ggplot2-animations-to-demonstrate-several-parallel-numerical-experiments-in-a-single-layout
#All errors are mine.
#This part needs imagemagick installed and will save the .gif to a temp directory based on your
#Imagemagick options.
##library(animation);
##saveMovie({#
## for (i in sort(unique(GDP.giant$year))) {#
## print(qplot(lambda,data=subset(GDP.giant,year==i),geom="density",fill=region,main = i)+opts(aspect.ratio = 1)+scale_x_continuous(limits=c(0,4)))#
## }#
##}, interval = 0.9, moviename = "lambda dist over time.gif", width = 600, height = 600);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment