Modern Toolmaking

caretEnsemble Classification example

2013-03-16T20:34:00.002-07:00

Here's a quick demo of how to fit a binary classification model with caretEnsemble. Please note that I haven't spent as much time debugging caretEnsemble for classification models, so there's probably more bugs than my last post. Also note that multi class models are not yet supported.

Right now, this code fails for me if I try a model like a nnet or an SVM for stacking, so there's clearly bugs to fix.

The greedy model relies 100% on the gbm, which makes sense as the gbm has an AUC of 1 on the training set. The linear model uses all of the models, and achieves an AUC of .5. This is a little weird, as the gbm, rf, SVN, and knn all achieve an AUC of close to 1.0 on the training set, and I would have expected the linear model to focus on these predictions. I'm not sure if this is a bug, or a failure of my stacking model.

New package for ensembling R models

2013-03-13T07:36:00.001-07:00

I've written a new R package called caretEnsemble for creating ensembles of caret models in R. It currently works well for regression models, and I've written some preliminary support for binary classification models.

At this point, I've got 2 different algorithms for combining models:

1. Greedy stepwise ensembles (returns a weight for each model)
2. Stacks of caret models

(You can also manually specify weights for a greedy ensemble)

The greedy algorithm is based on the work of Caruana et al., 2004, and inspired by the medley package here on github. The stacking algorithm simply builds a second caret model on top of the existing models (using their predictions as input), and employs all of the flexibility of the caret package.

All the models in the ensemble must use the same training/test folds. Both algorithms use the out-of-sample predictions to find the weights and train the stack. Here's a brief script demonstrating how to use the package:

Please feel free to submit any comments here or on github. I'd also be happy to include any patches you feel like submitting. In particular, I could use some help writing support for multi-class models, writing more tests, and fixing bugs.

Time series cross-validation 5

2013-01-24T16:22:00.001-08:00

The caret package for R now supports time series cross-validation! (Look for version 5.15-052 in the news file). You can use the createTimeSlices function to do time-series cross-validation with a fixed window, as well as a growing window. This function generates a list of indexes for the training set, as well as a list of indexes for the test set, which you can then pass to the trainControl object.

Caret does not currently support univariate time series models (like arima, auto.arima and ets), but perhaps that functionality is coming in the future? I'd also love to see someone write a timeSeriesSummary function for caret that calculates error at each horizon in the test set and a createTimeResamples function, perhaps using the Maximum Entropy Bootstrap.

Here's a quick demo of how you might use this new functionality:

Error metrics for multi-class problems in R: beyond Accuracy and Kappa

2012-07-06T10:08:00.001-07:00

The caret package for R provides a variety of error metrics for regression models and 2-class classification models, but only calculates Accuracy and Kappa for multi-class models. Therefore, I wrote the following function to allow caret:::train to calculate a wide variety of error metrics for multi-class problems:

This function was prompted by a question on cross-validated, asking what the optimal value of k is for a knn model fit to the iris dataset. I wanted to look at statistics besides accuracy and kappa, so I wrote a wrapper function for caret:::confusionMatrix and auc and logLoss from the Metric packages. Use the following code to fit a knn model to the iris dataset, aggregate all of the metrics, and save a plot for each metric to a pdf file:

This demonstrates that, depending on what metric you use, you will end up with a different model. For example, Accuracy seems to peak around 17:

While AUC and logLoss seem to peak around 6:

You can also increase the number of cross-validation repeats, or use a different method of re-sampling, such as bootstrap re-sampling.

Time series cross-validation 4: forecasting the S&P 500

2012-06-11T11:15:00.000-07:00

I finally got around to publishing my time series cross-validation package to github, and I plan to push it out to CRAN shortly.

You can clone the repo using github for mac, for windows, or linux, and then run the following script to check it out:
This script downloads monthly data for S&P 500 (adjusted for splits and dividends), and, for each month form 1995 to the present, fits a naive model, an auto.arima() model, and an ets() model to the past 5 year's worth of data and uses those models to predict S&P 500 prices for the next 12 months (note that the progress bar doesn't update if you register a parallel backend. I can't figure out how to fix this bug):

The naive model outperforms the arima and exponential smoothing models, both of which take into account seasonal patterns, trends, and mean-reversion! Furthermore, we're not just using any arima/exponential smoothing model: at each step we're selecting the best model, based on the last 5 years worth of data. (The ets model slightly outperforms the naive model at the 3 month horizon, but not the 2 month or 4 month horizons).

Forecasting equities prices is hard!

My first R package: parallel differential evolution

2012-01-23T06:53:00.000-08:00

UPDATE: a better parallel algorythm will be included in a future version of DEoptim, so I've removed my package from CRAN. You can still use the code from this post, but keep Josh's comments in mind.

Last night I was working on a difficult optimization problems, using the wonderful DEoptim package for R. Unfortunately, the optimization was taking a long time, so I thought I'd speed it up using a foreach loop, which resulted in the following function:

Here's what's going on: I divide the bounds for each parameter into n segments, and use a foreach loop to run DEoptim on each segment, collect the results of the loop, and then return the optimization results for the segment with the lowest value of the objective function. Additionally, I defined a "parDEoptim" class to make it easier to combine the results during the foreach loop. All of the work is still being done by the DEoptim algorithm. All I've done is split up the problem into several chunks.

Here is an example, straight out of the DEoptim documentation:

In theory, on a 20-core machine, this should run a bit faster than the serial example. Note that you may need to set itermax for the parallel run at a higher value than (itermax for the serial run)/(number of segments), as you want to make sure the algorithm can find the minimum of each segment. Also note that, in this example, there are 20 segments on the interval c(-10,-10) to c(10,10), which means that 2 of the segments have boundaries at c(1,1), which is the global minimum of the function. The DEoptim algorithm has no trouble finding a solution at the boundary of the parameter space, which is why it's so easy to parallelize.

Rumor has it that the next version of DEoptim will include foreach parallelization, but if you can't wait until then, I rolled up the above function into an R package and posted it to CRAN. Let me know what you think!

Benchmarking time series models

2011-12-29T10:47:00.000-08:00

This is a quick post on the importance of benchmarking time-series forecasts. First we need to reload the functions from my last few posts on times-series cross-validation. (I copied the relevant code at the bottom of this post so you don't have to find it).

Next, we need to load data for the S&P 500. To simplify things, and allow us to explore seasonality effects, I'm going to load monthly data, back to 1980.

The object "Data" has monthly closing prices for the S&P 500 back until 1980. Next, we cross validate 3 time series forecasting models: auto.arima, from the forecast package, a mean forecast, that returns the mean value over the last year, and a naive forecast, which assumes the next value of the series will be equal to the present value. These last 2 forecasts serve as benchmarks, to help determine if auto.arima would be useful for forecasting the S&P 500. Also note that I'm using BIC as a criteria for selecting arima models, and I have trace on so you can see the results of the model selection process.

After the 3 models finish cross-validating, it is useful to plot their forecast errors at different horizons. As you can see, auto.arima performs much better than the mean model, but is constantly worse than the naive model. This illustrates the importance of benchmarking forecasts. If you can't constantly beat a naive forecast, there's no reason to waste processing power on a useless model.

Finally, here is all the code in one place. Note that you can parallelize the cv.ts function by loading your favorite foreach backend.

Time series cross-validation 3

2011-12-12T08:37:00.000-08:00

I've updated my time-series cross validation algorithm to fix some bugs and allow for a possible xreg term. This allows for cross-validation of multivariate models, so long as they are specified as a function with the following paramters: x (the series to model), xreg (independent variables, optional), newxreg (xregs for the forecast), and h (the number of periods to forecast). Note that h should equal the number of rows in the xreg matrix. Also note that you need to forecast the xreg object BEFORE forecasting your x object. For example, if you wish to forecast 12 months into the future, your xreg object should have 12 extra rows.

Here is the source code for the new function:

And here is an example, using a linear model with xregs:

I am particularly excited about this code because it will allow me to apply arbitrary machine learning algorithms to forecasting problems. For example, I could create an xreg matrix of lags and use a support vector machine, neural network, or random forest to make 1-step forecasts. I am planning to release this code as a package on CRAN, once I finish the documentation. I'm also planning to re-work the function a bit to return an S3 class, containing:

Predicted values at each forecast horizon, including beyond the length of the input time series.
Actual values at each forecast horizon, for easy comparison to #1.
Matrix of average error at each horizon.
The final model.
Forecasts using the final model, from the last observation of x to the "max horizon".
A print method that will show #3.
A plot method that will plot #3.

Let me know if you have any suggestions or spot any bugs!

A pure R poker hand evaluator

2011-12-05T06:16:00.001-08:00

There's already a lot of great posts out there about poker hand evaluators, so I'll keep this short. Kenneth J. Shackleton recently released a very slick 5-card and 7-card poker hand evaluator called SpecialK. This evaluator is licensed under GPL 3, and is described in detail in 2 blog posts: part 1 and part 2. Since the provided code is open source, I felt free to hack around with it a bit, and ported the python source to R.

You can download my code from github, and save it to ~/SpecialK/R. Run the following script to initialize the evaluator and test it out. Higher numbers=better hands.

My R code is almost identical to the python source, as all I did was change python classes to R lists and make the lists 1-indexed (R) rather than 0-indexed (python). I also used the wonderful bitops package for bitwise operations. Obviously this code should be vectorized, but I don't have time to do that right now. I also used the "compiler" package to speed things up somewhat, but the "SevenEval" file still takes a long time to load.

Right now, I'm getting about 18,000 hands per second on my core i5 laptop, which is pretty crappy compared to the 250 million hands per second Shackleton reports for the C++ version. Vectorization should further increase performance, but I don't think pure R is ever going to approach C++ in terms of raw speed. While this was a fun port to make, I think this code is an obvious candidate for a re-write using the Rcpp package.

I'm interested to see how far the pure R code can be optimized.

Time series cross-validation 2

2011-11-22T07:12:00.001-08:00

In my previous post, I shared a function for parallel time-series cross-validation, based on Rob Hyndman's code. I thought I'd expand on that example a little bit, and share some additional wrapper functions I wrote to test other forecasting algorithms. Before you try this at home, be sure to load the cv.ts and tsSummary functions from my last post.

These functions add random walk models, the theta method, structural time series, seasonal decomposition, and simple mean forecasts to our cross-validation repertoire. The following code fits each of these models to the example dataset, and charts their accuracies out to a forecast horizon of 12 months. Note that none of this code runs in parallel, but if you wish, you can parallelize things by loading your favorite foreach backend. I would suggest running meanf, rwf, thetaf, and the linear model before loading a parallel backend, as all of these methods run very fast and do not need parallelization.

Here is the resulting figure. As you can see, the mean forecast is very inaccurate, but provides a useful baseline. The random walk forecast and the theta forecast are both improvements, but they ignore the function's seasonal component and have a clear seasonal error pattern. StructTS and STL are clustered down at the bottom with accuracies similar to the linear model, arima model, and exponential smoothing model:

If we ignore the mean, random walk, and theta forecasts, we get the following figure:

As you can see, the structural time series is close to the exponential smoothing model in accuracy, while the to seasonal decomposition models are consistently worse. The arima model still outperforms all other models, at every forecast horizon.

(Note that I've found a couple of bugs in my ts.cv function. It seems to not be working when fixed=TRUE, and it also doesn't like being told to just look at 1-step forecasts. I'll try to fix both bugs soon.)

Functional and Parallel time series cross-validation

2011-11-21T08:24:00.001-08:00

Rob Hyndman has a great post on his blog with example on how to cross-validate a time series model. The basic concept is simple: You start with a minimum number of observations (k), and fit a model (e.g. an arima model) to those observations. You then forecast out to a certain horizon (h), and compare your forecasts to the actual values for that series. You then add the next observation (k+1), and repeat the process until you run out of data. This gives you a matrix of forecast accuracies at various horizons (1 step ahead, 2 steps ahead, all the way to h steps ahead). You then take the mean of each column of the matrix, and get the model's average accuracy at that horizon. This method is analogous to leave-one-out cross-validation.

There are 2 variations to this method:
1. Use a fixed training window. In this case, when you add an observation to your "training" series, you drop the first observation, keeping the training window fixed.
2. Increment by n at each step, rather than 1. This is analogous to k-fold cross-validation. In this case, your forecast error is more unstable, and it's a good idea to average error across ALL horizons when evaluating the model.

This technique is very useful, because it allows you to define a horizon of interest (say 1 month or 12 months), and then asses how well your model performs at that horizon. Furthermore, you can use this data to compare various models, including different types of models, such as linear models vs. arima models vs. exponential smoothing model.

However, time series cross-validation is very time consuming, particularly for arima and exponential smoothing models. Therefore, I thought it would be a good idea to parallelize Hyndman's algorithm, using the foreach package in R. Furthermore, I wrapped the entire thing into a single function, which allows you to easily change the type of cross validation by altering the minObs (k), stepSize (n), and fixed-length or growing window parameters. My function takes an argument tsControl, which contains each of these parameters, as well a summary function to calculate your error metric (such as MAE). I've structured it similarly to the caret packages's train function.

The FUN argument should be a function that takes 2 parameters: x and h. x is a univariate time series, and h is a forecast horizon. The function should build a model using x, and then return h forecasts. Here are some examples of this function, for linear models, arima models, and exponential smoothing models. As you can see, they all return the same output: a vector of h point forecasts.

Here is my replication of Hyndman's example. Note that I create a "tsControl" list which contains the parameters for the cross-validation algorithm. This example is not parallelized, but you can easily change that by loading your favorite "foreach" backend, such as doMC or doRedis. Note that parallelization introduces a lot of overhead, and actually seems to slow down the linear model. I recomend only running the more complicated models (such as ets and Arima) in parallel.

This function produces the same results as Hyndman's example:

This code is very early stage, and could contain bugs. Please comment if you find one, or if you think of any ways to improve this algorithm. Feel free to try it out with your own forecasting functions!

Backtesting Part 4: random strategies

2011-10-21T07:21:00.000-07:00

Note: This post is NOT financial advice! This is just a fun way to explore some of the capabilities R has for importing and manipulating data.

In part 2, we found that our 200-day high, hold 100 days strategy yielded average annual returns on the S&P 500 index of 7% in a backtest going back to 1950. That's pretty good, but it's entirely possible that we got lucky and this was due to chance. One way to test this is to compare our strategy to a strategy that randomly chooses long or short each day. I tested just such a strategy on the S&P 500 index, again going back to 1950. I also adjusted the returns series for splits and dividends, and assumed trade costs were 0.5%. These are the exact same conditions I tested the "high-and-hold" strategy under.

Here is a histogram of average annual returns for these random strategies:

As you can see these strategies perform pretty dismally: the 0.5% trading cost quickly eats up all the available capital.

It would probably be a good idea to model a set dollar amount, and have trades costs a fixed amount, something like $7. It might also be a good idea to make the random strategies auto-correlated, which would probably better reflect actual investor behavior. The code for this test is below the fold, if anyone wants to modify it and make these improvements.

Backtesting a Simple Stock Trading Strategy: Part 3

2011-10-17T08:01:00.000-07:00

Note: This post is NOT financial advice! This is just a fun way to explore some of the capabilities R has for importing and manipulating data.

In a previous post, I examined a simple stock trading strategy: Find the high point over the last 200 days, and buy the stock if it's been less than 100 days since that high. Otherwise, have no position.

What if we use different parameters than 200-day high and hold 100 days? How will that affect our strategy? First of all, we have to reload the data for the S&P 500 index and re-define the functions used to implement our strategy.

Next, we must decide the range of parameters we wish the test for our strategy. I've decided to use a "grid search" to thoroughly examine the parameter space. Somewhat arbitrarily, I've decided to test the values from 5-500, by 5, for both parameters. This gives us 100 possible values for each parameter, or 10000 total. Good thing the "daysSinceHigh" function is pretty fast!

Because my processing power is limited, I'm only going to look at every 5th value in this parameter space. The first order of business is to calculate a matrix containing each n-Day high series, where the first column is the number of days since the 5-day high, the second column is the number of days since the 10-day high, etc. This matrix has 100 columns:

Next, I make a list with 100 elements. Each element represents a holding period, which I will apply to a copy of the "n-Day high matrix" from the previous step. For example, the 1st element in the list is a matrix representing a 5-day holding period. The first column in this matrix represents buying at the 5-day high, and holding for 5 days. This is equivalent to buy-and-hold. The second column represents buying at the 10-day high, and holding for 5 days. The third column represents buying at the 15-day high and so on. I repeat this process for each element in the 100-matrix list, which gives us an object representing every possible permutation of our strategy.

It is then a relatively easy thing to calculate the returns associated with each permutation of the strategy, by using the "sweep" function to multiply each column of each matrix by the daily returns for our stock

Now we have a list of matrices of returns. Each column of a matrix represents the returns of our strategy, using a different set of parameters. This allows us to calculate cumulative returns for each set of parameters, and make a nifty graph that shows the relationship between nHigh, nHold, and returns.

This graph uses a custom color ramp function, which was created by Andrie on StackOverflow. The color of each point in the corresponds to how high the returns are at that point. The X axis is number of days to use for the nHigh, and the yAxis is the number of days to use for nHold. As you can see, 100 days seems to be a solid holding period across many values of nHigh, but by using a different value of nHigh, we could increase returns substantially.

Of course, just because these values worked in the past doesn't mean they will work in the future. Still, it's good to see that our arbitrary parameters (which performed well in the last post), fall inside a wide range of parameters that yield a positive return for our strategy. This brings up an interesting question: how DO we select parameters for our strategy? How can we tell how well our parameter selection strategy would have performed in the past, given that we've optimized our selection based on of our knowledge of the past?

For homework, think about how overfitting and cross-validatation apply to this problem...

BONUS CODE: This creates some nifty 3D charts, using the rgl library.

Recession forecasting III: A Better Naive Forecast

2011-09-20T11:42:00.000-07:00

In Recession Forecasting Part II, I compared the accuracy of Hussman's recession forecasts to the accuracy of a naive forecast that assumed the current state of the recession variable would continue next month.

An anonymous commentator pointed out a great error in the naive forecast, which is that the NBER BACKDATES their recession indicator. In other words, a recession may start, but you would not know about it for up to 6 months. Therefore, assuming the naive forecast would know of the recession in the second month is akin to predicting the future.

This creates the need for a naive forecast that only has access to the information available at a given point in time. Fortunarly, our anonymous commentator kindly provided this data, which he constructed using the recession announcement dates from the NBER. Previous to 1980, he assumed a 6-month lag in identifying turning points. This data is available as an excel spreadsheet (with commentary) and a csv file (suitable for import into R).

Using this modified dataset, we can construct truly naive recession forecast, which assumes that the current known state of the recession variable will continue indefinitely.

This naive forecast is 66.9% accurate, which is pretty good. This number also make a lot more intuitive sense than the 97.6% accuracy of my previous naive forecast. Compared to this new metric, Hussman's accuracy of 81.6% is actually quiet good. However, keep in mind that Hussman's "future" knowledge of past recessions may have biased his selection of variables and parameters.

Finally, Hussman's model uses Stock Prices, PMI, the yield curve, credit spreads, and employment as inputs. I was wondering if there are any significant lags in any of these variables that could introduce look-ahead bias into my implementation of Hussman's model.

Backtesting Part 2: Splits, Dividends, Trading Costs and Log Plots

2011-09-16T11:03:00.000-07:00

Note: This post is NOT financial advice! This is just a fun way to explore some of the capabilities R has for importing and manipulating data.

UPDATE 9/21/2011: Costas correctly points out that I should lag my strategy vector by one day, as today's returns are determined the position we chose yesterday. I updated the code, results, and charts to reflect this. Please alert me to any similar errors in the future.

In my last post, I demonstrated how to backtest a simple momentum-based stock trading strategy in R. However, there were a few issues with my implementation: I ignored splits, dividends and transaction costs, all of which can have a large impact on a strategy. I also came up with a better plot to help show how a given strategy compares to a benchmark, and wrapped everything together into one function.

First of all, we need to reload our functions from the last post. These functions define our strategy and analyze its performance.

Next, we can test our strategy. I've added a couple new indexes:

This function does a lot of lifting:
1. It loads the data and adjusts the closing price for splits and dividends. It uses the splits/dividends data from yahoo, but performs its own, more accurate calculations.
2. It determines a position series, based on the "daysSinceHigh" function. This part is the same as in my last post.
3. It determines trades, which are defined as days when today's position is different from the previous day's positions. I assumed that transactions costs are 0.5% of equity, so on trading days I subtracted 0.005 from my Returns.
4. It makes a plot. This plot is different from the charts.PerformanceSummary we used last time. The first plot shows cumulative returns of my strategy and the index, while the second plot shows the relative performance of my strategy over the benchmark (also known as alpha). The third plot shows drawdowns.
5. It returns a data table of statistics, comparing the strategy to the benchmark.

I tested this strategy on GSPC, FTSE, DJI, N225, EEM, EFA, and GLD. (The last 3 are ETFs). The strategy performs well on some indexes, and poorly on others. Here's the results of my backtest:

As you can see, this strategy tends to reduce drawdowns, but it also sometimes reduces overall returns. In some cases, you could leverage up the strategy, which would increase both returns and drawdowns, but that's the subject of another post.

Here's a buncha charts:

Correlations among US Stocks: Is it really time to fire your adviser?

2011-09-15T14:35:00.000-07:00

Note: This post is NOT financial advice! This is just a fun way to explore some of the capabilities R has for importing and manipulating data.

The Financial Times says it's time to "Fire your Adviser" because correlations among US stocks are at their highest levels since the financial crisis. Unfortunately, they only provide data going back 3 months and it's in a boring table, rather than an awesome chart. After reading this article, I immediately pulled out quantmod and PerfomanceAnalytics in R. I used quantmod to download the daily data (from yahoo finance) for each index listed in the article, and then used PerformanceAnalytics to analyze and graph it. After and hour or so of fiddling around, I ended up with this chart:

Each line on this chart represents the rolling correlation between a major asset class and the S&P 500 (represented by SPY). I used a 90-day window to calculate the correlations, so each point on each line is looking backwards at a 90-day timeframe. As you can see, correlations have indeed gone up in the last few months, but there have been other periods in 2010 and 2011 with similarly high correlations.

This next chart looks at just the major US sectors, and takes the average of each of their correlations with SPY. This shows the overall correlation of the US stock market since 2007. Again, things are highly correlated right now, but we've been here before.

Here is the code I used to generate these charts. Feel free to comment on my implementation, and I'll be happy to make any improvements and update the charts.

Backtesting a Simple Stock Trading Strategy

2011-09-13T06:06:00.000-07:00

Note: This post is NOT financial advice! This is just a fun way to explore some of the capabilities R has for importing and manipulating data.

UPDATE 9/21/2011: Costas correctly points out that I should lag my strategy vector by one day, as today's returns are determined the position we chose yesterday. I updated the code, results, and charts to reflect this. Please alert me to any similar errors in the future.

I recently read a post on ETF Prophet that explored an interesting stock trading strategy in Excel. The strategy is simple: Find the high point of the stock over the last 200 days, and count the number of days that have elapsed since that high. If its been more less than 100 days, own the stock. If it's been more than 100 days, don't own it. This strategy is very simple, but it yields some impressive results. (Note; however, that this example uses data that has not been adjusted from splits or dividends and could contain other errors. Furthermore, we're ignoring trading costs and execution delays, both of which affect strategy performance.)

Implementing this strategy in R is simple, and provides numerous advantages over excel, the primary of which is that pulling stock market data into R is easy, and we can test this strategy on a wide range of indexes with relatively little effort.

First of all, we download data for GSPC using quantmod. (GSPC stands for the S&P 500 index). Next, we construct a function to calculate the number of days since the n-day high in a time series, and a function to implement our trading strategy. The latter function takes 2 parameters: the n-day high you wish to use, and the numbers of days past that high you will hold the stock. The example is 200 and 100, but you could easily change this to the 500-day high and see what happens if you hold the stock 300 days past the high before bailing out. Since this function is parameterized, we can easily test many other versions of our strategy. We pad the beginning of our strategy with zeros so it will be the same length as our input data. (If you wish for a more detailed explaination of the daysSinceHigh function, see the discussion on cross-validated).

We multiply our position (0,1) vector by the returns from the index to get our strategy's returns. Now we construct a function to return some statistics about a trading strategy, and compare our strategy to the benchmark. Somewhat arbitrarily, I've decided to look at cumulative return, mean annual return, sharpe ratio, winning %, mean annual volatility, max drawdown, and max length drawdown. Other stats would be easy to implement.

Results:

As you can see, this strategy compares favorably to the default "buy-and-hold" approach.

Finally, we test our strategy on 3 other indexes: FTSE which represents Ireland and the UK, the Dow Jones Industrial Index, which goes back to 1896, and the N225, which represents Japan. I've functionalized the entire process, so you can test each new strategy with 1 line of code:

Results:

FTSE:

DJIA:

N225:

The strategy out-performs the other indexes as well. It even performs better than the N225 index, mainly by staying out of it.

Feel free to test this strategy with other parameters, or on other indexes. For homework, think of possible ways that I have fooled myself in this backtest, and post them in the comments. One example of this is that we have't looked at transaction costs, which might be significant...

25+ more ways to bring data into R

2011-08-26T14:07:00.000-07:00

The rdatamarket post on the Revolutions blog and this post on Decision Stats reminded me about my list of Data APIs/feeds available as packages in R on Cross-Validated (which is a great site that you all should use). Many of these packages are from Omega Hat, which is an awesome site.

This is the most comprehensive list I'm aware of, so please alert me to any ommisions:

Free Data:
Data Source - Package
Google Finance historical data - quantmod
Google Finance balance sheets - quantmod
Yahoo Finance historical data - quantmod
Yahoo Finance historical data - tseries
Yahoo Finance current options chain - quantmod
Yahoo Finance historical analyst estimates - fImport
Yahoo Finance current key stats - fImport - seems to be broken
CFE historic futures prices - qmao
OANDA historic exchange rates/metal prices - quantmod
FRED historic macroeconomic indicators - quantmod
World Bank historic macroeconomic indicators - WDI
Google Trends historic search volume data - RGoogleTrends
Google Docs - RGoogleDocs
Twitter - twitteR
Zillow - Zillow
New York Times - RNYTimes
US Census 2000 - UScensus2000
infochimps - infochimps
datamarket - rdatamarket - requires free account
Factual.com - factualR
Geocode addresses - RDSTK
Map coordinates to political boundaries - RDSTK
Weather Underground - Roll your own
Google News - Roll your own
Earth Sciences netCDF Data - Roll your own
Climate Data - Roll your own
Public health data - Roll your own
FishBase - rfishbase

Paid Data:
Bloomberg - RBloomberg
LIM - LIM
Trades and Quotes from NYSE - RTAQ
Interactive Brokers - IBrokers

Useful Tools:
RCurl
RJSON
RJSONIO
XML
scraper
digitizer

Because it's Friday: Spurious correlation edition

2011-08-26T08:26:00.000-07:00

If the Flight of the Concords taught me anything, it's that you can't trust Australians. This morning I was poking around the DataMarket site, when I noticed something suspicious about Australian sheep production:

I decided to investigate further: just what are the Australians doing with all those poor sheep? How might this nefarious plot impact the rest of the world? I was shocked when I noticed an apparent relationship with Global Warming:

Hot on the trail of something big, I pulled both datasets into R, using the rdatamarket package, and ran a quick correlation, which confirmed my worst suspicions:

The Great Australian Sheep Decline is clearly the cause of Global Warming! To quantify this effect, I built a simple linear model, which proved to be highly signifigant:

For every 100 Million Sheep that disappear from Australia, the Earth warms by 0.71 degrees C! Given that there are still 85 million sheep left in Australia, the potential for further warming is gigantic!

=========================================================================================
In all seriousness, Australian sheep have about as much to do with Global Warming as pirates. Sheep numbers have been declining in New Zealand, Urugay, Iceland, and many other countries as well. If you wish to explore other spurious correlations, here's the R code I used for this post:

Graphically analyzing variable interactions in R

2011-08-23T10:30:00.000-07:00

I studied Ecology as an undergraduate, which meant I spent a lot of time gathering and analyzing field data. One of the basic tools we used to look for relationships in a large set of variables was correlation and scatterplot matrices. Each of these requires a single line of code in R:

The 'pairs' function in R contains a lot of additional options, which can be used to make very informative plots. These options can get a little cumbersome, but fortunately several package authors have written wrapper functions that automatically enable some extra magic. Two such packages are psych and PerformanceAnalytics. I happen to prefer the 1 liner from PerformanceAnalytics, but it's a matter of personal taste:

This chart contains a LOT of information: On the diagonal are the univariate distributions, plotted as histograms and kernel density plots. On the right of the diagonal are the pair-wise correlations, with red stars signifying significance levels. As the correlations get bigger the font size of the coefficient gets bigger. On the left side of the diagonal is the scatter-plot matrix, with loess smoothers in red to help illustrate the underlying relationship. This is one of my favorite plots in R, because it combines a large amount of information into one command and one easy to follow plot. In fact, this plot contains more information than is revealed by the 1st two commands in this post!

Of course, you can use this command on data from other domains besides Ecology. PerformanceAnalytics is intended for the analysis of financial data, so lets put it through its paces. First we download some financial data (a stock index, a bond index, and a gold index) from yahoo finance using quantmod, and then combine the daily close series of those indexes into one dataframe. I'm not 100% happy with the legend in the plot, but I wanted to show how the correlations between these indexes have changed over the years. I also skipped red (color #2) in the plots and in the legend, because the loess smoother is also red.

Finally, I'd like to acknowledge Stephen Turner over on cross-validated for inspiring this post.

Recession forecasting II: Assessing Hussman's Accuracy

2011-08-22T09:01:00.000-07:00

In my last post on recessions, I implemented John Hussman's Recession Warning Composite in R. In this post I will examine how well this index performs and discuss how we might improve it. If you would like to follow along at home, be sure to run the code from the last post, before running anything from this post.

First of all, lets evaluate how predictive Hussman's index is of recessions next month:

This code simply compares the current value of USREC (US Recessions) to last month's value of the recession warning composite. By this measure, the recession warning composite is only 81.55% accurate, with a 95% confidence interval of [75%,87%].

Next, let's evaluate a warning ANYTIME in the last 6 months to the current value of USREC.

By this measure, the forecast is even worse: the accuracy is 73.62% [66%,80%]. Interestingly, this measure has a very high 'Negative Predictive Value' (.9896), which indicates if the recession warning composite has been 0 for the past 6 months, you can be reasonable sure there will be no recession this month.

Finally, let's make a naive recession forecast, and predict that the current value of USREC will be equal to it's previous value:

This forecast is 97.62% accurate! [94%,99%]. Therefore, I have to conclude that Hussman's recession warning composite, while interesting to implement, is not particularly useful for forecasting recessions. However, it may be that Hussman is primarily concerned with forecasting when recessions START and END. Given that the current state of US recessions is highly predictive of the next state of US recessions, this might be a valid approach. Still, I'm struggling to find a useful way of interpreting Hussman's index.

If you have any ideas for using or improving Hussman's index, please leave a comment.

Using the google prediction API from R

2011-08-10T11:16:00.000-07:00

Google has a "black box" prediction API that they provide for use with creating recommender systems or filtering spam. Furthermore, they provide an R package for interfacing that API, but try as I might I cannot get it to work under windows. Here are the instructions for setting up the API to run in R under linux. I haven't tried this out yet, so let me know in the comments if it works, or if you can get it to run on Windows.

First we have to setup the Google Prediction API, as well as some dependencies:
1. Go to the Google APIs Console. This is your home base for managing google APIs.
2. In the upper left hand corner of the website (under the Google APIs logo) is a dropdown menu. Use this to create a new project, called something informative like "R predictions."
3. Activate the Google storage API and turn it on. Activating may require opening a new page.
4. Activate the Google prediction API and turn it on. Activating may require opening a new page.
5. Click on the "Billing" tab, and make sure billing is enabled. You may have to enter your billing information. Note that you get 5GB of free storage through the end of 2011, and there's a free quota on the prediction API for 5MB trained per day and 100 predictions per day, up to 20,000 total predictions.
6. Click the "Google Storage" tab, and make a note of the "x-goog-project-id." You will need this when installing GSUtils.

Next we have to install some software on our computer to enable communication between R and the prediction API:
1. Install python, if you do not already have it.
2. Make sure you can run python from the command prompt. You may need to add python to your "path" or "environment" variables to do this. On windows, run the command prompt as administrator.
3. Install the R packages rjson and RCurl using install.packages() in R.
4. Make sure you can open .tar archives. This is no problem on Mac/Linux systems, but on windows you need 7zip.
5. Download GSUtil, and follow the directions to install it on your system. This is the tricky part.
6. When you run GSUtil for the first time, make sure to use the following command: python gsutil config -b to allow gsutil to open a web page and authorize access to your google storage account.
6. When prompted, enter the project ID you recorded in part 1.
7. Download the googlepredictionapi package.
8. Open R, and setwd() to the folder containing the downloaded package.
9. Install the R package from source using this command:

Now we're all set to start using the prediction API:
1. First we need to create a bucket to store our data. Do this from the Google Storage Web Console. Name your bucket something useful, like rdata. Don't use capital letters or symbols.
2. Run the following script to test that everything works. Note that you have to save your data frame as a .csv file before GSUtil can upload it to google storage for modeling:

Good luck! Here are some links for future reference:
1. Google directions for installing the googlepredictionapi package in R
2. Google directions for installing gsutil
3. Google API console for managing APIs and billing
4. Google storage console for managing buckets
5. Google APIs overview/introduction

Scraping web data in R

2011-08-10T06:36:00.000-07:00

In my last post, I went through a lot of effort to scrape the PMI index off the ISM website. It turns out that was unnecessary effort, as commentator "senne" pointed out that this index is available from FRED, with the symbol NAPM. I've updated my code, which now pulls all the data straight from FRED.

However, it was surprisingly easy to scrape web data into R, using the readHTMLTable function in the XML package. I thought I'd keep the code I used on my blog, as it's a good example of how easily you can pull web data into R.

Forecasting recessions

2011-08-09T09:29:00.000-07:00

John Hussman has a Recession Warning Composite that I am attempting to replicate/improve. The underlying data seems to be easy enough to get from FRED using the quantmod package in R. I don't quite understand the index Hussman is using for commercial paper, so I used the '3-month AA financial commercial paper index' from FRED.

The PMI index requires scraping a table from the ISM website, which is easy enough to do with the XML package.

Here's my code so far, please leave a comment and let me know what you think:

Parallel random forests using foreach

2011-07-22T06:56:00.000-07:00

There's been some discussion on the kaggle forums and on a few blogs about various ways to parallelize random forests, so I thought I'd add my thoughts on the issue.

Here's my version of the 'parRF' function, which is based on the elegant version in the foreach vignette:

This function works very simply: you pass it a vector of mtry values, and it fits a random forest using each of those values and returns the combined result. You can all pass any additional parameters you want (like ntree) to the randomForest function.

I think this functions provides 2 improvements over previous implementations. #1 is you can use any parallel backend you want. doRedis is my current favorite, as it's cross-platform and fault-tolerant and let's me commandeer idle laptops around the house/office when a random forest is taking too long to fit. #2 is the argument .inorder=FALSE in the foreach function, which provides a small performance improvement as it lets R combine the random forests as they finish, rather than forcing R to combine them in the order they start.

Lets say you want a random forest with 5000 trees. The default value for ntree is 500, so we use rep(4,10) as the argument for the function.

Maybe we're unsure of the optimal mtry value, and want combine 2 ensembles of 2500 trees. Then we use the argument c(rep(3,5),rep(4,5)). This gives us 2500 trees with mtry=3 and 2500 with mtry=4. I like to think of this as a sort of meta-ensemble of decision trees, but I've yet to see it improve my predictive accuracy.

At the very least, this can help with those damn 'out of memory' errors I've been getting on my laptop when fitting random forests to large datasets.