Blog Archives

Downloading SP500 stock price data with Julia

[Update]

In order to provide a more convenient and stable way of getting data, I meanwhile started the Julia package EconDatasets.jl. It basically bundles many download routines for different datasets in a common interface, and is more easy to keep maintain through automated testing on travis. With this package, you also can download the same data with just two lines of code:

Pkg.add("EconDatasets")
using EconDatasets
# download data into package subdirectory data 
getDataset("SP500")
# load data into session workspace
data = dataset("SP500")

Original post

In a previous post I already described how to download stock price data for all constituents of the SP500 with R. Meanwhile, however, I shifted all of my coding to the new and very promising technical computing language julia. I obviously wanted to also enjoy the convenience of downloading stock price data directly in julia. Hence, in this post I will describe how one could use julia to download stock price data from Yahoo!Finance in general, and – as an advanced use case – how to download adjusted closing prices for all SP500 constituents. We also will make use of julia’s parallel computing capabilities, in order to tremendously speed up the downloading and processing steps compared to my original code in R.

Before we start, we first need to make sure that we have all relevant packages installed. As a format for storing time series data I always rely on the TimeData package from the official package repository:

Pkg.add("TimeData")

The functions required for data download are bundled in the not yet registered EconDatasets package, which we can clone from github:

Pkg.clone("https://github.com/JuliaFinMetriX/EconDatasets.jl.git")

1 Download data using readYahooFinance

Now that all required packages are installed, we will first demonstrate a simple use case for the low-level function readYahooFinance. This functions allows easy access to the complete stock data provided by Yahoo!Finance, comprising fields open, high, low, close, volume and adjusted close. As an example, we now download data for NASDAQ-100, S&P 500 and EURO STOXX 50 Index.

using TimeData
using Dates
using EconDatasets

tickerSymbs = ["^NDX"
               "^GSPC"
               "^STOXX50E"]


dates = Date(1960,1,1):Date(2014,7,20)

indexData = [readYahooFinance(dates, symb) for symb in tickerSymbs]

## display first five dates of NASDAQ
indexData[1][1:5, :]
idx Open High Low Close Volume Adj_Close
1985-10-01 221.24 224.32 221.13 224.28 153160000 112.14
1985-10-02 224.28 225.08 221.56 221.65 164640000 110.82
1985-10-03 221.68 222.37 220.24 221.74 147300000 110.87
1985-10-04 221.74 221.74 219.71 220.15 147900000 110.07
1985-10-07 220.15 220.27 216.35 216.4 128640000 108.2

Using comprehension, the function readYahooFinance can be applied to all ticker symbols in an Array. The output will be an Array of type Any, with individual entries being of type Timematr.

When we focus on one variable for each stock only, we can store the data more concisely in a single TimeData object. Therefore, we join individual stocks at the their idx entries. We do not want to lose any data at this step, so that we will use an outer join in order to get a row for each date that occurs for at least one of the individual stocks. Missing values will be replaced by NA, so that we now get an object of type Timenum, as Timematr objects are not allowed to contain NAs.

adjCloseData = indexData[1][:Adj_Close]
names!(adjCloseData.vals, [symbol(tickerSymbs[1])])

for ii=2:3
    nextStock = indexData[ii] |>
         x -> x[:Adj_Close]
    names!(nextStock.vals, [symbol(tickerSymbs[ii])])
    adjCloseData = joinSortedIdx_outer(adjCloseData, nextStock)
end

adjCloseData[[1:10; (end-10):end], :]
idx ^NDX ^GSPC ^STOXX50E
1950-01-03 NA 16.66 NA
1950-01-04 NA 16.85 NA
1950-01-05 NA 16.93 NA
1950-01-06 NA 16.98 NA
1950-01-09 NA 17.08 NA
1950-01-10 NA 17.03 NA
1950-01-11 NA 17.09 NA
1950-01-12 NA 16.76 NA
1950-01-13 NA 16.67 NA
1950-01-16 NA 16.72 NA
2014-07-04 NA NA 3270.47
2014-07-07 3910.71 1977.65 3230.92
2014-07-08 3864.07 1963.71 3184.38
2014-07-09 3892.91 1972.83 3203.1
2014-07-10 3880.04 1964.68 3150.59
2014-07-11 3904.58 1967.57 3157.05
2014-07-14 3929.46 1977.1 3185.86
2014-07-15 3914.46 1973.28 3153.75
2014-07-16 3932.33 1981.57 3202.94
2014-07-17 3878.01 1958.12 3157.82
2014-07-18 3939.89 1978.22 3164.21

2 Download adjusted closing prices of SP500 constituents

Now that we already have seen a first use case of function readYahooFinance, we now want to try the capabilities of julia and the EconDatasets package with a more challenging task. Hence, we want to download adjusted stock prices for all constituents of the SP500 in a fully automated way. Therefore, we first need to get a list of the ticker symbols of all constituents, which we can get from the S&P homepage. However, this list is stored as an Excel sheet with .xls extension, and we need to read in this binary file with package Taro.

To make Taro work, you first need to make sure that it is able to find Java on your system. If your path deviates from the default settings, just make sure to set the respective JAVA_LIB environment variable in your .bashrc file. In my case, the variable is set as follows:

# 64-bit machine
# JAVA_LIB="/usr/lib/jvm/java-7-openjdk-amd64/jre/lib/amd64/server/"

# 32-bit machine
JAVA_LIB="/usr/lib/jvm/java-7-openjdk-i386/jre/lib/i386/server/"
export JAVA_LIB

We can now install and load package Taro:

Pkg.add("Taro")
using Taro
Taro.init()
Found libjvm @ /usr/lib/jvm/java-7-openjdk-i386/jre/lib/i386/server/

If something with your Taro configuration is not correct, you will get an error at this step. In this case, you could simply download and export the Excel sheet to .csv manually, which you then can read in with function readtable from the DataFrames package.

Otherwise, you can use Taro to download and read in the respective part of the Excel sheet:

url = "http://us.spindices.com/idsexport/file.xls?hostIdentifier=48190c8c-42c4-46af-8d1a-0cd5db894797&selectedModule=Constituents&selectedSubModule=ConstituentsFullList&indexId=340"
filepath = download(url)
constituents = Taro.readxl(filepath, "Constituents", "A10:B511")
head(constituents)
Constituent Symbol
3M Co MMM
Abbott Laboratories ABT
AbbVie Inc. ABBV
Accenture plc ACN
ACE Limited ACE
Actavis plc ACT

We now should have name and ticker symbol of each SP500 constituent stored as a DataFrame. In my case, however, there even is one ticker symbol too much, although I do not know why:

(nTicker, nVars) = size(constituents)
501
2

An inconsistency that I will not further invest at this point. In addition, however, some of the ticker symbols are automatically read in as boolean values, and we will have to convert them to strings first. Let’s display all constituents with boolean values:

isBoolTicker = [isa(tickerSymbol, Bool) for tickerSymbol in
                constituents[:Symbol]]

constituents[find(isBoolTicker), :]
Constituent Symbol
AT&T Inc true
Ford Motor Co false

The reason for this is that the respective ticker symbols are “T” and “F”, which will be interpreted as boolean values. Once we did correct for this mistake, we transform the array of ticker symbols into an Array of type ASCIIString.

indTrue = find(constituents[2] .== true)
indFalse = find(constituents[2] .== false)

constituents[indTrue, 2] = "T"
constituents[indFalse, 2] = "F"

tickerSymb = ASCIIString[constituents[:Symbol]...]
tickerSymb[1:5]
MMM
ABT
ABBV
ACN
ACE

Now that we already have a list of all ticker symbols, in principle we could apply the same procedure as before: download each stock, extract the adjusted closing prices, and join all individual price series. However, as we have 500 stocks, this procedure would already take approximately 15 minutes if each individual stock took only 2 seconds. Hence, we strive for a much faster result using julia’s parallel computing capabilities, and this is already implemented as function readYahooAdjClose.

Under the hood, readYahooAdjClose uses a map-reduce structure. As the map step, for any given ticker symbol we download the data, extract the adjusted closing prices and rename the column to its ticker symbol. As reduce step we need to specify some operation that combines the individual results of the map step – in our case, this is function joinSortedIdx_outer.

Let’s now set the stage for parallel computation, add three additional processes and load the required packages on each process.

addprocs(3)

@everywhere using Dates
@everywhere using DataFrames
@everywhere using TimeData
@everywhere using EconDatasets
3-element Array{Any,1}:
 2
 3
 4

To run the parallelized code, simply call function readYahooAdjClose:

dates = Date(1960,1,1):Date(2014,7,20)

## measure time
t0 = time()

@time vals = readYahooAdjClose(dates, tickerSymb, :d)

t1 = time()
elapsedTime = t1-t0
mins, secs = divrem(elapsedTime, 60)

Downloading of all 500 stocks did only take:

println("elapsed time: ", int(mins), " minutes, ", ceil(secs), " seconds")
elapsed time: 3 minutes, 45.0 seconds

Now we convert the data of type Timedata to type Timenum and store the result in the EconDatasets data directory:

valsTn = convert(Timenum, vals)
pathToStore = joinpath(Pkg.dir("EconDatasets"), "data", "SP500.csv")
writeTimedata(pathToStore, valsTn)

3 Visualize missing values

From previous experiences I already know that the saying “you get what you pay for” also holds for the free Yahoo!Finance database: the data comes with a lot of missing values. In order to get a feeling for the data quality, we want to visualize all missing values. Therefore, we sort all assets with respect to the number of missing values.

(nObs, nStocks) = size(valsTn)

NAorNot = Array(Int, nObs, nStocks)

for ii=1:nStocks
    NAorNot[:, ii] = isna(vals.vals[ii])*1
end

nNAs = sum(NAorNot, 1)
p = sortperm(nNAs[:])
tickSorted = tickerSymb[p]
NAorNotSorted = NAorNot[:, p]

For the visualization, we rely on plot.ly, as we get an interactive graphic this way. This allows the identification of individual stocks in the graphic by simple mouse clicks. If you want to replicate the figure, you will need to sign in with an own and free plot.ly account.

Pkg.clone("https://github.com/plotly/Plotly.jl")
using Plotly
## sign in with your account
## Plotly.signin("username", "authentication")

nToPlot = nStocks
datsStrings = ASCIIString[string(dat, " 00:00:00") for dat in valsTn.idx]
data = [
        [
         "z" => NAorNotSorted[:, 1:nToPlot],
         "y" => tickSorted[1:nToPlot],
         "x" => datsStrings, 
         "type" => "heatmap"
         ]
        ]

response = Plotly.plot([data], ["filename" => "Published graphics/SP500 missing values", "fileopt" => "overwrite"])
plot_url = response["url"]

Although the resulting graphic easily could be directly embedded into this html file, you will need to follow this link to watch it, since the plot.ly graphic is a large data file and hence takes quite some time to load. Nevertheless, I also exported a .png version of the graphic, which you can find below. Thereby, red dots are representing missing values.

missing_values.png

4 Get logarithmic returns

So now we already have our closing price data at our hands. However, in financial econometrics we usually analyze return data, as prices are non-stationary almost always. Hence, we now want to derive returns from our price series.

For this step I rely on the also not yet registered Econometrics package:

Pkg.clone("https://github.com/JuliaFinMetriX/Econometrics.jl.git")
using Econometrics

The reason for this is that I can make use of function price2ret then, which implements a slightly more sophisticated approach to return calculation. This can best be described through an example, which we can set up using function testcase from the TimeData package:

logPrices = testcase(Timenum, 4)
idx prices1 prices2
2010-01-01 100 110
2010-01-02 120 120
2010-01-03 140 NA
2010-01-04 170 130
2010-01-05 200 150

As you can see, the example logarithmic prices have a single missing observation at January 3rd. Straightforward application of the logarithmic return calculation formula hence will result in two missing values in the return series:

logRetSimple = logPrices[2:end, :] .- logPrices[1:(end-1),:]
idx prices1 prices2
2010-01-02 20 10
2010-01-03 20 NA
2010-01-04 30 NA
2010-01-05 30 20

In contrast, function price2ret assumes that single missing values in the middle of a price series are truly non-existent: there is no observation, because the stock exchange was closed and there simply was no trading. For example, this could easily happen if you have stocks from multiple countries with different holidays. Note, that this kind of missingness is different to a case where the stock exchange was open and trading did occur, but we were not able to observe the resulting price (for a more elaborate discussion on this point take a look at my blog post about missing stock price data). Using function price2ret, you ultimately will end up with a stock return series with only one NA for this example:

## get real prices
logRetEconometrics = price2ret(logPrices; log = true)
idx prices1 prices2
2010-01-02 20 10
2010-01-03 20 NA
2010-01-04 30 10
2010-01-05 30 20

So the final step in our logarithmic return calculation is to apply price2ret to the logarithmic prices, specify the usage of differences for return calculations through log = true, and multiply the results by 100 in order to get percentage returns.

percentLogRet = price2ret(log(valsTn); log = true).*100

Alternatively, you could get discrete returns with:

percentDiscRet = price2ret(valsTn; log = false).*100
Advertisements

Coping with missing stock price data

1 Missing stock price data

When downloading historic stock price data it happens quite frequently that some observations in the middle of the sample are missing. Hence the question: how should we cope with this? There are several ways how we could process the data, each approach coming with its own advantages and disadvantages, and we want to compare some of the most common approaches in this post.

In any case, however, we want to treat missing values as NA and not as Julia’s built-in NaN (a short justification on why NA is more suitable can be found here). Hence, data best should be treated through DataFrames or – if the data comes with time information – through type Timenum from the TimeData package. In the following, we will use these packages in order to show some common approaches to deal with missing stock price data using an artificially made up data set that shall represent logarthmic prices.

The reason why we are looking at logarthmic prices and returns instead of normal prices and net returns is just that logarithmic returns are defined as simple difference between logarithmic prices of successive days:

\displaystyle r_{t}^{\log}=\log(P_{t}) - \log(P_{t-1})

This way, our calculations simply involve nicer numbers, and all results equally hold for normal prices and returns as well.

We will be using the following exemplary data set of logarthmic prices for the comparison of different approaches:

using TimeData
using Dates
using Econometrics

df = DataFrame()
df[:stock1] = @data([100, 120, 140, 170, 200])
df[:stock2] = @data([110, 120, NA, 130, 150])

dats = [Date(2010, 1, 1):Date(2010, 1, 5)]

originalPrices = Timenum(df, dats)
idx stock1 stock2
2010-01-01 100 110
2010-01-02 120 120
2010-01-03 140 NA
2010-01-04 170 130
2010-01-05 200 150

One possible explanation for such a pattern in the data could be that the two stocks are from different countries, and only the country of the second stock has a holiday at January the 3rd.

Quite often in such a situation, people just refrain from any deviations from the basic calculation formula of logarithmic returns and calculate the associated returns as simple differences. This way, however, each missing observation NA will lead to two NAs in the return series:

simpleDiffRets = originalPrices[2:end, :] .- originalPrices[1:(end-1), :]
idx stock1 stock2
2010-01-02 20 10
2010-01-03 20 NA
2010-01-04 30 NA
2010-01-05 30 20

For example, this also is the approach followed by the PerformanceAnalytics package in R:

library(tseries)
library(PerformanceAnalytics)

stockPrices1 <- c(100, 120, 140, 170, 200)
stockPrices2 <- c(110, 120, NA, 130, 150)

## combine in matrix and name columns and rows
stockPrices <- cbind(stockPrices1, stockPrices2)
dates <- seq(as.Date('2010-1-1'),by='days',length=5)
colnames(stockPrices) <- c("A", "B")
rownames(stockPrices) <- as.character(dates)
(stockPrices)

returns = Return.calculate(exp(stockPrices), method="compound")
nil nil
20 10
20 nil
30 nil
30 20

When we calculate returns as the difference between successive closing prices P_{t} and P_{t-1}, a single return simply represents all price movements that happened at day t, including the opening auction that determines the very first price at this day.

Thinking about returns this way, it obviously makes sense to assign a value of NA to each day of the return series where a stock exchange was closed due to holiday, since there simply are no stock price movements on that day. But why would we set the next day’s return to NA as well?

In other words, we should distinguish between two different cases of NA values for our prices:

  1. NA occurs because the stock exchange was closed this day and hence there never were any price movements
  2. the stock exchange was open that day, and in reality there were some price changes. However, due to a deficiency of our data set, we do not know the price of the respective day.

For the second case, we really would like to have two consecutive values of NA in our return series. Knowing only the prices in t and t+2, there is no way how we could reasonably deduce the value that the price did take on in t+1. Hence, there are infinitely many possibilities to allocate a certain two-day price increase to two returns.

For the first case, however, it seems unnecessarily rigid to force the return series to have two NA values: allocating all of the two-day price increase to the one day where the stock exchange was open, and a missing value NA to the day that the stock exchange was closed doesn’t seem to be too arbitrary.

This is how returns are calculated by default in the (not yet registered) Econometrics package.

logRet = price2ret(originalPrices, log = true)
idx stock1 stock2
2010-01-02 20 10
2010-01-03 20 NA
2010-01-04 30 10
2010-01-05 30 20

And, the other way round, aggregating the return series again will also keep NAs for the respective days, but otherwise perform the desired aggregation. Without specified initial prices, aggregated prices will all start with value 0 for logarithmic prices, and hence express something like normalized prices that allow a nice comparison of different stock price evolutions.

normedPrices = ret2price(logRet, log = true)
idx stock1 stock2
2010-01-01 0 0
2010-01-02 20 10
2010-01-03 40 NA
2010-01-04 70 20
2010-01-05 100 40

To regain the complete price series (together with a definitely correct starting date), one simply needs to additionally specify the original starting prices.

truePrices = ret2price(logRet, originalPrices, log = true)
idx stock1 stock2
2010-01-01 100 110
2010-01-02 120 120
2010-01-03 140 NA
2010-01-04 170 130
2010-01-05 200 150

In some cases, however, missing values NA may not be allowed. This could be a likelihood function that requires real values only, or some plotting function. For these cases NAs easily can be removed through imputation. For log price plots, a meaningful way would be:

impute!(truePrices, "last")
idx stock1 stock2
2010-01-01 100 110
2010-01-02 120 120
2010-01-03 140 120
2010-01-04 170 130
2010-01-05 200 150

However, for log returns, the associated transformation then would artificially incorporate values of 0:

impute!(logRet, "zero")
idx stock1 stock2
2010-01-02 20 10
2010-01-03 20 0
2010-01-04 30 10
2010-01-05 30 20

As an alternative to replacing NA values, one could also simply remove the respective dates from the data set. Again, there are two options how this could be done.

First, one could remove any missing observations directly in the price series:

originalPrices2 = deepcopy(originalPrices)
noNAprices = convert(TimeData.Timematr, originalPrices2, true)
idx stock1 stock2
2010-01-01 100 110
2010-01-02 120 120
2010-01-04 170 130
2010-01-05 200 150

For the return series, however, we then have a – probably large – multi-day price jump that seems to be a single-day return. In our example, we suddenly observe a return of 50 for the first stock.

logRet = price2ret(noNAprices)
idx stock1 stock2
2010-01-02 20 10
2010-01-04 50 10
2010-01-05 30 20

A second way to eliminate NAs would be to remove them from the return series.

logRet = price2ret(originalPrices)
noNAlogRet = convert(TimeData.Timematr, logRet, true)
idx stock1 stock2
2010-01-02 20 10
2010-01-04 30 10
2010-01-05 30 20

However, deriving the associated price series for this processed return series will then lead to deviating end prices:

noNAprices = ret2price(noNAlogRet, originalPrices)
idx stock1 stock2
2010-01-01 100 110
2010-01-02 120 120
2010-01-04 150 130
2010-01-05 180 150

as opposed to the real end prices

originalPrices
idx stock1 stock2
2010-01-01 100 110
2010-01-02 120 120
2010-01-03 140 NA
2010-01-04 170 130
2010-01-05 200 150

The first stock now suddenly ends with a price of only 180 instead of 200.

2 Summary

The first step when facing a missing price observation is to think whether it might make sense to treat only one return as missing, assigning the complete price movement to the other return. This is perfectly reasonable for days where the stock market really was closed. In all other cases, however, it might make more sense to calculate logarithmic returns as simple differences, leading to two NAs in the return series.

Once there are NA values present, we can chose between three options.

2.1 Keeping NAs

Keeping NA values might be cumbersome in some situations, since some functions might only be working for data without NA values.

2.2 Replacing NAs

In cases where NAs may not be present, there sometimes might exist ways of replacing them that perfectly make sense. However, manually replacing observations in some way means messing with the original data, and one should be cautious to not incorporate any artificial patterns this way.

2.3 Removing NAs

Obviously, when dates with NA values for only some variables are eliminated completely, we simultaneously lose data for those variables where observations originally were present. Furthermore, eliminating cases with NAs for returns will lead to price evolutions that are different to the original prices.

2.4 Overview

Possible prices:

idx simple diffs single NA replace \w 0 rm NA price rm NA return
2010-01-01 100, 110 100, 110 100, 100 100, 110 100, 110
2010-01-02 120, 120 120, 120 120, 120 120, 120 120, 120
2010-01-03 140, NA 140, NA 140, 120    
2010-01-04 170, 130 170, 130 170, 130 170, 130 150, 130
2010-01-05 200, 150 200, 150 200, 150 200, 150 180, 150

Possible returns:

idx simple diffs single NA replace \w 0 rm NA price rm NA return
2010-01-02 20, 10 20, 10 20, 10 20, 10 20, 10
2010-01-03 20, NA 20, NA 20, 0    
2010-01-04 30, NA 30, 10 30, 10 50, 10 30, 10
2010-01-05 30, 20 30, 20 30, 20 30, 20 30, 20

Downloading SP500 stock price data with R

1 Summary

In this blog post I want to show how to download data of all constituents of the S&P500 index, which is one of the leading stock market indices for US equity, from Yahoo!Finance. According to Standard and Poor’s, “the index includes 500 leading companies and captures approximately 80% coverage of available market capitalization.” The target will be to get a clean and broadly applicable dataset, which then can be used for many types of financial applications, such as risk management and asset management.

Instead of running the R code snippets on your own, you can also directly download the most recent version of the resulting dataset which I keep for my own research: historic prices and clean return series. These datasets, however, could slightly deviate from the results obtained in this blog post, since they probably build on a more recent version of the input files.

Update: Due to the Yahoo!Finance data disclaimer, I will not provide the application-ready and cleaned-up final dataset. I am sorry, but in order to get the data, you need to run the steps on your own.

2nd Update: Meanwhile I did replicate the whole exercise in my new statistical software language of choice: julia. Not only does this dramatically speed up the procedure, but I also did process the data with regards to missing values in a slightly improved manner. You can check out this new version here, and in addition you can find a comparison of different approaches to deal with missing values at this blog post.

2 Data considerations

The best way to test the adequacy of a model is testing its performance on out-of-sample data. With this application in mind, the dataset should be large enough so that it can be divided into a part that is used to fit the model and another part that can be used for out-of-sample testing. The problem with financial data is that asset return distributions vary over time, thereby posing additional challenges to risk and asset management. This most easily can be seen in times of crisis, where markets clearly exhibit higher levels of volatility. Any econometric model now should be able to replicate reality at any given point in time, and therefore needs to be able to capture both calm and turbulent market periods. Therefore, your model ideally should also be tested on both market states, so that your out-of-sample period should entail a crisis period as well. Similarly, of course, your model ideally should be fitted to both types of market states as well. For example, fitting a model to calm market times only will most likely lead to improper results during financial crisis, since crisis would be one region of the unconditional distribution that the model has never seen before. For a very extreme example, think about fitting a model to positive returns only. You would clearly fail to capture the overall unconditional distribution, and not to speak of the conditional distribution of negative returns.

According to this argumentation, we will pick January the 1st, 1995, as the beginning of the data sample. The data sample hence entails two financial crisis: the bursting of the dot-com bubble in 2001, and the 2008 financial crisis following the bursting of the US housing bubble. The first crisis now could be used to fit some econometric model, while we can test its performance in both risk or asset management applications out-of-sample during the second crisis.

3 Data download

The first step – of course – is downloading the data from Yahoo!Finance. This means that we need to get all the individual assets that comprise the index. A list of all current constituents can be accessed on the S&P homepage, and can be directly downloaded as xls-file here, or – if the direct link is broken – by manually clicking on the full constituents list on the S&P 500 homepage. However, to make things easier, you can download this already processed csv-file, which then can be directly read into R. Following my convention, you should copy the file in the subfolder raw_data of the project’s directory, to indicate that this list comprises externally created data.

In the past, this list was where I did encounter the first problem with the data: the list of 500 components was already missing the symbol of one stock. A problem that seems to be fixed by now.

After reading in the list of constituents, the function get.hist.quote from the tseries package is used to download the associated data in a for-loop. Note, that I did not conduct any linearization or preallocation with this code, so that the code’s processing time could be further decreased in general. Anyways, since individual historic price series need to be merged to common dates, acceleration most likely would result in significantly less readable code.

## downloads historic prices for all constituents of SP500
library(zoo)
library(tseries)                        

## read in list of constituents, with company name in first column and
## ticker symbol in second column
spComp <- read.csv("raw_data/sp500_constituents.csv" ) 

## specify time period
dateStart <- "1995-01-01"               
dateEnd <- "2013-05-01"

## extract symbols and number of iterations
symbols <- spComp[, 2]
nAss <- length(symbols)

## download data on first stock as zoo object
z <- get.hist.quote(instrument = symbols[1], start = dateStart,
                    end = dateEnd, quote = "AdjClose",
                    retclass = "zoo", quiet = T)

## use ticker symbol as column name 
dimnames(z)[[2]] <- as.character(symbols[1])

## download remaining assets in for loop
for (i in 2:nAss) {
   ## display progress by showing the current iteration step
   cat("Downloading ", i, " out of ", nAss , "\n")

   result <- try(x <- get.hist.quote(instrument = symbols[i],
                                     start = dateStart,
                                     end = dateEnd, quote = "AdjClose",
                                     retclass = "zoo", quiet = T))
   if(class(result) == "try-error") {
      next
   }
   else {
      dimnames(x)[[2]] <- as.character(symbols[i])

      ## merge with already downloaded data to get assets on same dates 
      z <- merge(z, x)                      

   }


}

## save data
write.zoo(z, file = "data/all_sp500_price_data.csv", index.name = "time")

Although this code did work for me, please note that it is not implemented in a completely robust way, as the download of the first asset is not executed in a try block. Hence, if problems occur with this asset, you might encounter an error.

4 Visualizing missing values

The best way to identify any data errors with such high-dimensional data always is through visualization. As first step, we will identify assets and dates with missing values. One easy way of looking at this is to order assets with respect to the number of occurring NAs, and highlighting entries with missing values. Usually, I tend to create my visualizations with ggplot2 whenever possible. However, given the size of the dataset, this will be quite time-consuming. Hence, I will also implement a comparable version of the visualization using only low-level graphical tools of R with faster execution. Missing values are shown in light blue and pink respectively.

## load libraries
library(zoo)
library(ggplot2)
library(lattice)
library(reshape)

## read in sp500 prices
z <- read.zoo(file = "data/all_sp500_price_data.csv", header = T)
nAss <- ncol(z)

## decode missing values with numeric ones
miss <- is.na(coredata(z))*1            # multiplication converts
                                        # logical to numeric matrix

## summing up missing values per column and ordering with respect to
## number of missing values 
ranks <- rank(colSums(miss), ties.method = "random")
miss[, ranks] <- miss

################################
## visualization with ggplot2 ##
################################

## missing values as dataframe 
missDf <- data.frame(miss)

## append time index as first column
missDf$time <- factor(index(z))
missDf <- missDf[c(nAss+1, 1:nAss)]

## processing dataframe for ggplot2
missMt <- melt(missDf, id.vars = c("time"))

## plotting values
p <- (ggplot(missMt, aes(x = time, y = variable, colour = value)) +
      geom_tile(aes(fill = value)))

## getting rid of date factor labels
p2 <- p + theme(axis.line = element_blank(), 
                axis.text.x = element_blank(),
                axis.text.y = element_blank(), 
                axis.ticks = element_blank(),
                legend.position = "none",
                panel.background = element_blank(),
                panel.border = element_blank(), 
                panel.grid.major = element_blank(),
                panel.grid.minor = element_blank(),
                plot.background = element_blank()) +
  xlab("time") + ylab("stocks")

##################################
## visualization with levelplot ##
##################################

## getting rid of column names
colnames(miss) <- NULL
levelplot(miss, aspect = 1, col.regions = cm.colors,
          xlab = "time", ylab = "Stocks", main = "Missing values",
          colorkey = F)

missing_values-1.png

missing_values-2

As can be seen, there are a significant number of missing values for about one fifth of the 500 S&P components. In addition, there are also a few dates where data is missing for a substantial share of stocks. We now want to simply eliminate stocks with multiple missing values, since we assume that these data problems are only arising due to deficiencies of the data provider. In other words, it generally should be possible to get clean data on the chosen sample period through usage of a proprietary and more comprehensive database. We hence do not tackle stocks with large ratios of missing values in any ad-hoc way, since we do not want to include any spurious effects into the subsequent analysis.

As threshold value we choose to exclude all assets with more than nine missing values, and afterwards any dates with more than two missing values. Eliminated dates are shown together with the number of missing observations in a Google Charts html table. After completely removing stocks and dates with multiple NAs, all occasionally occurring remaining NAs will simply be set to 0 in the associated return series. For example, for a case of one missing observation for one asset, like

incompleteData <- cbind(c(4, NA, 6), c(6, 9, NA))
4 6
nil 9
6 nil

this procedure will lead to a return of zero for the first pair of days, and a return corresponding to the full price change for the second pair of days.

missingValues <- which(is.na(incompleteData))

filledValues <- incompleteData
filledValues[missingValues] <- incompleteData[missingValues-1] 
(filledValues)
(diff(filledValues))
     [,1] [,2]
[1,]    4    6
[2,]    4    9
[3,]    6    9
     [,1] [,2]
[1,]    0    3
[2,]    2    0

At the end, we will have a matrix of approximately 4600 observations for each of the approximately 350 remaining stocks.

## load packages
library(zoo)
library(googleVis)
library(fGarch)

rm(list=ls())

## read sp500 prices
z <- read.zoo(file = "data/all_sp500_price_data.csv", header = T)

######################################################
## eliminate stocks with more than 2 missing values ##
######################################################

## get indices of stocks with max 9 missing values
miss <- is.na(coredata(z))*1
indices <- colSums(miss) < 10           
suitableStocks <- z[, indices]

## show number of remaining stocks
print(length(suitableStocks[1, ]))

## create table with missing values to show dates with many NAs
miss <- is.na(coredata(suitableStocks))*1
datesWithNAs <- rowSums(miss)
missingValues <- data.frame(dates = index(z)[datesWithNAs > 3],
                            number = datesWithNAs[datesWithNAs > 3])

## show missing values in html file
p <- gvisTable(missingValues)
write(p$html$chart, file = "pics/table_of_missingValues.html")

####################################################
## remove dates with more than two missing values ##
####################################################

datesToKeep <- rowSums(miss) < 3
pricesNaDatesAndStocksRemoved <- suitableStocks[datesToKeep, ]


#############################################################
## calculate percentage log returns and replace NAs with 0 ##
#############################################################

## show number of remaining NAs
print(sum(is.na(pricesNaDatesAndStocksRemoved)))

## avoid logical indexing with zoo object
allPrices <- coredata(pricesNaDatesAndStocksRemoved)
missingValues <- which(is.na(allPrices))

## fill missing values with value of previous day
filledValues <- allPrices
filledValues[missingValues] <- allPrices[missingValues - 1]
logRetNumeric <- 100*diff(log(filledValues))

## create zoo object again
returnDates <- index(pricesNaDatesAndStocksRemoved)[-1]
logRet <- zoo(x = logRetNumeric, order.by = returnDates)

## check that no NAs exist anymore
if(sum(is.na(logRet)) != 0){
   stop("still NAs in return series!!")
   ## could be due to NAs on consecutive days
}

## round to five decimal places
logRet <- round(logRet, 5)

## store data
write.zoo(logRet, file = "data/all_sp500_clean_logRet.csv")

Again, at this point the code is not implemented in a completely robust manner. If you have an NA at the first day of the sample, you would try to access an element at index 0. This, of course, would lead to an error. Furthermore, I do not deal explicitly with the case of consecutive NAs. If this case occurred, you would be warned, as the last if-statement would throw an exception. For both cases one should rather easily find some go-around with for loops. However, this most likely would kill your performance. And any more sophisticated solution I did skip in order to keep the code simple.

5 Zoo object pitfalls

When removing all of the remaining NAs, bear in mind that zoo objects can not be accessed by logical indexing! Hence, the following assignment does not work:

logRet[is.na(logRet)] <- 0

Instead, we need:

coredata(logRet)[is.na(logRet)] <- 0