Exploratory Data Analysis

When encountering a data set for the first time, it is often a good idea to explore the distribution of each variable.

Here is a simple R function that will create a nice display to aid in exploring the univariate distribution of a numeric vector.

For example, using the USArrests dataset:

data(USArrests)
eda.plots(USArrests$Murder, "Murder arrests (per 100,000)")

Source code for the eda.plots function:

#--------------------------------------------------------------------
# Univariate Exploratory Data Analysis
#--------------------------------------------------------------------
# PANELS
#--------------------------------------------------------------------
# (1) Top Left: Histogram
# (2) Top Right: Boxplot
# (3) Bottom Left: Density Plot
# (4) Bottom Right: Normal QQ-Plot
#--------------------------------------------------------------------
# ARGUMENT  | INPUT            | DESCRIPTION
#--------------------------------------------------------------------
# x         | numeric vector   | variable to be analyzed
# vname     | character vector | variable label & filename
# outdir    | character vector | output directory
# dev       | character vector | graphics device
# overwrite | logical vector   | overwrite existing files or not
# col       | character vector | color of graph elements
#--------------------------------------------------------------------
eda.plots <- function(x,
					  vname="edaplots",
					  outdir="./",
					  dev="png",
					  overwrite=TRUE,
					  col="black") {
	
	##  check inputs
	
	if(is.vector(x) && class(x)!="numeric") {
		stop("x is not a numeric vector\n")
	} else {
		stop("x is not a vector\n")
	}
	
	char.vars <- c(vname,outdir,dev,col)
	for(v in char.vars) {
		if(class(v)=="character") {
			if(nchar(v)>50)
				stop(sprintf("%s is too long",v))
			if(length(v)>1) {
				v <- v[1]
				warning(sprintf("%s has been coerced to its first element",v))
			}
		} else {
			stop(sprintf("%s is not a character vector",v))
		}
	}
	
	if(class(overwrite)=="logical") {
		if(length(overwrite)>1) {
			overwrite <- overwrite[1]
			warning(sprintf("%s has been coerced to its first element",v))
		}
	} else {
		stop("overwrite is not a logical vector")
	}
	
	## create directory if needed
	
	dir.create(file.path(outdir), showWarnings=FALSE)
	fname <- gsub("[[:blank:][:punct:]]", "", vname)
	plot.name <- sprintf("%s/%s.%s", outdir, tolower(fname), tolower(dev))
	
	if(!overwrite) {
		if(any(grepl(plot.name,dir())))
			stop(sprintf("%s currently exists already and will not be overwritten",
						 plot.name))
	}

	## plots
	
	# set up layout
	do.call(tolower(dev), list(plot.name))
	
	par(mfrow=c(2,2), mai=c(0.75, 0.75, 0.25, 0.25))
					# mai=c(bottom,left,top,right)
	
	# graph parameters
	cex.axis <- 1.25
	cex.lab  <- 1.5
	cex      <- 1.5
	lwd      <- 2
	
	# histogram - top left panel
	h.col <- ifelse(col=="black", "white", col)
	hist(x,
		 col=h.col,
		 lwd=lwd,
		 main="",
		 xlab="",
		 cex.axis=cex.axis,
		 cex.lab=cex.lab
	)
	
	# boxplot - top right panel
	b.col <- ifelse(col=="black", "white", col)
	boxplot(x,
			horizontal=TRUE,
			lwd=lwd,
			col=b.col,
			ylab=vname,
			cex.axis=cex.axis,
			cex.lab=cex.lab
	)
	
	# densityplot - bottom left panel
	plot(density(x),
		 col=col,
		 lwd=lwd,
		 main="",
		 xlab="",
		 cex.axis=cex.axis,
		 cex.lab=cex.lab
	)
	rug(x, col=col)
	
	# QQ - bottom right panel
	qqnorm(x,
		   col=col,
		   cex=cex,
		   lwd=lwd,
		   main="",
		   cex.axis=cex.axis,
		   cex.lab=cex.lab
	)
	qqline(x)
	
	dev.off()
	invisible()
}

Case Study: State of the Union Twitter Sentiment

During the 2013 State of the Union address, I thought it would be an interesting experiment to collect tweets containing the reserved hashtag ‘#SOTU’.

I wrote an R script that collected the latest 1000 tweets containing ‘#SOTU’ every 30 seconds. While this program was running, I took notes on the times and issues President Barack Obama spoke on. Using these two datasets, I was able to put together this plot of average tweet sentiment over time with the volume of tweets below.

From this information, we can see which issues had positive and negative reactions.

Positive Reaction:

• Health Care
• Education
• Immigration
• Minimum Wage
• Cyber Threat

Negative Reaction:

• Green Technology
• Afghanistan
• Gun Control

To learn to query data and quantify sentiment, please check out my other posts:
Extracting Data From Twitter
Scoring Sentiment of Tweets

Twitter Sentiment Analysis

Following my last post Extracting Data from Twitter, we are now ready to to play with the tweets we grabbed from Twitter.

The first thing we will do is authenticate our connection with Twitter.

load("twitCred.RData")
registerTwitterOAuth(twitCred)

Now, we need a way to score how ‘good’ or ‘bad’ each tweet is. This is actually an interesting problem. I will be using an adapted solution developed from Jeffrey Breen.

# load files containing positive & negative words
positive.words = scan("positive-words.txt",
					  what='character',
					  comment.char=';')
negative.words = scan("negative-words.txt",
					  what='character',
					  comment.char=';')

# dependent function
tryTolower = function(x){

	# create missing value
	# this is where the returned value will be
	y = NA

	# tryCatch error
	try_error = tryCatch(tolower(x), error=function(e) e)

	# if not an error
	if (!inherits(try_error, "error"))
	y = tolower(x)

	return(y)
}

# sentiment score function
score.sentiment = function(sentences, pos.words, neg.words, .progress='none'){

    # We got a vector of sentences.
	# plyr will handle a list or a vector as an "l" for us.
    # We want a simple array ("a") of scores back, so we use
    # "l" + "a" + "ply" = "laply":
    scores = laply(sentences, function(sentence, pos.words, neg.words){

        # Clean Up Sentences With R's Regex-driven Global Substitute, gsub():
        sentence = gsub('[[:punct:]]', '', sentence)
        sentence = gsub('[[:cntrl:]]', '', sentence)
        sentence = gsub('\\d+', '', sentence)

        # Convert to Lower Case:
		sentence = iconv(sentence, 'UTF-8', 'ASCII')
		sapply(sentence, function(x) tryTolower(x))

        # Split into Words. str_split is in the stringr package
        word.list = str_split(sentence, '\\s+')

        # sometimes a list() is one level of hierarchy too much
        words = unlist(word.list)

        # Compare Our Words to the Dictionaries of Positive & Negative Terms
        pos.matches = match(words, pos.words)
        neg.matches = match(words, neg.words)

        # match() returns the position of the matched term or NA
        # we just want a TRUE/FALSE:
        pos.matches = !is.na(pos.matches)
        neg.matches = !is.na(neg.matches)

        # and conveniently enough, TRUE/FALSE will be treated as 1/0 by sum():
        score = sum(pos.matches) - sum(neg.matches)

        return(score)
    }, pos.words, neg.words, .progress=.progress )

	scores.df = data.frame(score=scores, text=sentences)
	return(scores.df)
}

# function to grab tweets
grab.tweet <- function(search="#rstats", n=10) {

	# check search parameter
	if (!is.character(search)) {
		warning(sprintf("Search argument %s is not character.",search))
	}

	# check n parameter
	if (!is.numeric(n)) {
		warning(sprintf("Number %s is not numeric.",n))
	}

	# get tweets
	tweets <- searchTwitter(search, n=n, cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))

	# put in data frame
	tweets <- twListToDF(tweets)

	# output to temporary file and read back in to weed out erroneous characters
	write.csv(tweets, file="tmp_tweet.csv", row.names=FALSE)
	tweets <- read.csv(file="tmp_tweet.csv")

	# score tweets
	tweet.scores <- score.sentiment(tweets$text, positive.words, negative.words)

	# merge back to original data
	tweet.scores <- merge(tweets,tweet.scores)
	return(tweet.scores)
}

Now we can load tweets and analyze the sentiment score.

# get 500 latest tweets containing "#heat"
tweets <- grab.tweet(search="#heat", n=500)

# examine data frame with tweets and score
head(tweets)

# tabulate scores
table(tweets$score)

# plot scores and output as png file
require(lattice)
png(file="heat_bar.png", width=1000, height=750)
barchart(as.factor(tweets$score),
		 horizontal=FALSE,
		 col="red",
		 scales=list(x=list(cex=1.5), y=list(cex=1.5)),
		 xlab=list("Score",cex=2),
		 ylab=list("Frequency",cex=2),
		 main=list("Miami Heat Tweet Sentiment",cex=4))
dev.off()

This was run right after the Miami Heat won game 5 of the 2013 NBA Playoffs against the Indiana Pacers in the Eastern Conference finals.

Extracting Data from Twitter

Data mining has become an interesting hobby for me. The idea of collecting real-time data seems very powerful and useful. Using R, I will show you how to query tweets from Twitter.

Step 1:
Login in to Twitter Developers. Once you are in, create a new application. Click on create access token. Now click on settings and change the Application Type to Read, Write and Access Direct Messages and Allow this application to be used to sign on Twitter.

Step 2:
Start up R and run this script with your Consumer key and Consumer secret.

require("twitteR")

Key      <- "<< Consumer key >>"   # Consumer key
Secret   <- "<< Consumer secret >>" # Consumer secret
twitCred <- OAuthFactory$new(consumerKey=Key,
							 consumerSecret=Secret,
							 requestURL="https://api.twitter.com/oauth/request_token",
							 accessURL="https://api.twitter.com/oauth/access_token", 
							 authURL="https://api.twitter.com/oauth/authorize")
twitCred$handshake(cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))

R will spit out a URL to visit. Copy and paste this into your browser. Click authorize app and copy the PIN. Go back to your R console and paste the PIN and press enter.

Step 3:
Save credentials as .RData to bypass this step later.

save(twitCred,file="twitCred.RData")

This will save a file called “twitCred.RData” to your working directory.

Step 4:
Now every time you want to extract tweets from Twitter, all you need to do is load your credentials and authenticate.

load("twitCred.RData")
registerTwitterOAuth(twitCred)

Step 5:
Now we are ready to extract tweets from Twitter.

rstats <- searchTwitter('#rstats', n=100, cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))
rstatsDF <- twListToDF(rstats)

Markowitz Mean-Variance Portfolio

I have been meaning to re balance my portfolio for a while now and thought I would try a new investing strategy.  As a former math and economics student, I will pay tribute to the Nobel prize winning economist Harry Markowitz and try to construct a minimum variance portfolio using equities in the S&P 500.  To get the ticker symbols for each stock in the S&P 500, I went and downloaded the table from Wiki as a csv file.

Let’s begin by setting up our workspace:

# Set Up Workspace
rm(list=ls(all=TRUE))

# Load Required Packages
require("quantmod")
require("corpcor")
require("tawny")
require("quadprog")

# Read in all S&P 500 Companies
sp500 <- read.csv("../data/sp500.csv")

Now, query daily stock prices for each equity (this takes a while, go get some coffee):

# Input Stock Tickers to Run Through
tickers <- as.character(sp500$Ticker.symbol)

# Get Stock Data from and to a certain date
getSymbols(tickers, from="2013-01-01")

Calculate daily stock returns for each equity:

# Calculate Daily Returns and put in a data frame
DailyReturns <- do.call(merge, lapply(tickers, function(x){ periodReturn(get(x), period='daily') }))
names(DailyReturns) <- paste(tickers,".Return",sep="")

And solve the famous Markowitz problem. This will return the weights of our minimum variance portfolio:

# Assume returns are a matrix of log returns, then generate optimization inputs:
covar <- cov.shrink(DailyReturns)
N     <- ncol(DailyReturns)
zeros <- array(0, dim = c(N,1))

# Evaluate the optimization to generate minimum variance portfolio without short selling constraint
aMat <- t(array(1, dim = c(1,N)))
res  <- solve.QP(covar, zeros, aMat, bvec=1, meq = 1)

# Or, similar optimization with short selling constraint (i.e. non-negative weights)
aMat <- cbind(t(array(1, dim = c(1,N))), diag(N))
b0   <- as.matrix(c(1, rep.int(0,N)))
res  <- solve.QP(covar, zeros, aMat, bvec=b0, meq = 1)

# Finally, return portfolio attributes (similar to portfolio.optim)
y <- DailyReturns %*% res$solution
port <- list(pw = round(res$solution,3), px = y, pm = mean(y), ps = apply(y, 2, sd))

Finally, summarize this portfolio visually:

# Set Up Weights for Graphing
portfolio <- data.frame(as.character(tickers),port$pw)
names(portfolio) <- c("tick","wt")

# Remove Zeros
portfolio <- subset(portfolio, portfolio$wt!=0)

# Merge on Names
portfolio <- merge(portfolio, sp500, by.x="tick", by.y="Ticker.symbol", all.x=TRUE)
portfolio[4] <- NULL; portfolio[7] <- NULL

# Tabulate Sectors for Coloring
sectors <- data.frame(table(portfolio$GICS.Sector))
names(sectors) <- c("sector","count")

# Generate Colors
sectors <- cbind(sectors,rainbow(nrow(sectors)))
names(sectors)[3] <- "cols"

# Merge on colors
portfolio <- merge(portfolio, sectors, by.x="GICS.Sector", by.y="sector")

# Graph
png(file="../graphs/h-minvarpf.png", width=1500, height=1000)
par(mar=c(6, 6, 4, 2) + 0.1)

barplot(portfolio$wt,
		names=portfolio$tick, cex.names=1.75, las=2,
		main="Minimum Variance Portfolio Weights", cex.main=4,
		ylim=c(0,0.12), cex.axis=2, col=portfolio$cols)

legend("top", legend=unique(portfolio$GICS.Sector),
	   pch=15, col=unique(portfolio$cols), inset=0.025, cex=2, pt.cex=4)

dev.off()

Keep in mind that past performance is not indicative of future results.

Code was adapted from Quantitivity.