Pages

Sunday, March 3, 2013

What can we use the VIX for?


In part 1, we took a look at VIX and the relationship it had between historical volatility and realized volatility.

Continuing on, I thought I would take a look at next day returns and the VIX. There is a relationship between SPX and VIX in that when SPX drops, VIX typically rises. This leads the question: does a high VIX serve as a useful indicator of what might happen to SPX the next day?

A problem in answering this is quantifying what we mean by a "high" VIX. A VIX of 25 might be high if it has been under 20 for the last 3 months, but might be considered low if it has been around 30.

To deal with this I use a proportional measure based on its previous trading range. If it is at the top of its range, it gets a 1, and at the bottom it gets a 0. Anywhere in between it gets a number between 0 and 1.

I used data from 2000-2013 and a 252 period lookback, which is roughly 1 year, recording the next day return of SPX. I split this data into 5 groups, the first with relative VIX readings up to 0.2, the second between 0.2 and 0.4, and so on up to 1.0. Then I plotted the next day returns as a boxplot, giving the distribution of returns for a VIX reading in its relative quintile.



The box on the far left corresponds to the lowest relative VIX levels, while the box on the far right corresponds to the highest.

One thing that jumps out to me is that the mean of all these returns is more or less zero. Using this relative measure of VIX has no value as a predictor of next day returns, as whatever measure we use, a mean return of zero will see us end up flat (or in practice, probably down a bit). If there was a real edge, we would expect the mean return to be somewhere other than zero.

The second thing that I notice is the distribution of returns is a lot wider for the right most two boxes, which corresponds to relative VIX measure above 0.60, i.e. relatively higher volatility. Now this makes sense, as we would reasonably expect periods of higher volatility to have higher volatility.

It seems that high readings of VIX can serve as a useful predictor of higher periods of next day volatility, and for someone trying to implement a volatility filter, that could be quite useful.

Conclusion


To recap, we have seen that:

1) VIX is not very good at forecasting 30 day realized volatility.
2) Historical volatility can be a reasonable proxy for IV/VIX.
3) Higher periods of volatility do portend continued higher volatility, at least for the short term (1 day).

Note that 3 does contradict the mean reversion of volatility, because it is looking at the next day returns only. Over a longer period we would certainly expect it mean revert, otherwise we'd be experiencing some very, very bumpy rides ...

Also, I did similar analysis using RUT & RVX and NDX & VXN, both of which had very similar results. We will take a look at them in the next post where we use the findings from this post and the previous to implement a volatility filter.

Code for this post and the previous is up: here

Digging into the VIX


I wanted to revisit using some sort of volatility filter for systematic trading. In particular, if we are trading SPX, can we somehow use the VIX to produce better risk adjust returns? This is not about trading volatility, but more about using additional "out of band" data in our systems.

I thought I would take a look at the VIX, what is it, and can it help us?

In theory, the VIX is the market consensus of what future 30 day realized volatility of SPX will actually be, as derived from option prices of SPX. Anyone who watches or trades SPX has seen that a drop in SPX usually results a rise in VIX, or the IV in SPX option chains.

Occasionally this relationship breaks down, which can be a useful as an indicator in itself, but I had a sneaking suspicion that the VIX would have more in common with historical vol than what realized vol turns out to be over the next 30 calendar days. That is, the VIX is not particularly useful as a longer term forward looking indicator of volatility, at least no more than current historical volatility.

I did a bit of data gathering and calculated realized and historical volatility for SPX. I used 21 periods when calculating vols as the VIX covers 30 calendar days, which is roughly a month, and there are roughly 21 trading days in an average month.

First of all, we can take a look in hindsight at how well the VIX forecast of volatility matched the actual realized volatility over the next 21 days.



In this chart, the red line is what realized volatility turned out to be, going forward from a given day. It is impossible to know this in advance, so we are engaging in some serious data snooping, for educational purposes only. The teal line is the VIX. Below them is a plot of VIX - SPX realized vol, and flagrant abuse of gradients.

You can see that the VIX gives an estimate that is typically higher than what actually eventuated, a well-known phenomenon, which for convenience I will put in "The Volatility Risk Premium" basket.

There is one big divergence around the end of July, 2011, when the VIX was significantly underestimating what realized vol would actually turn out to be.  In August 2011 there was a large selloff and you can see the VIX jumped to a rather ebullient level.

If we moved the red line forward 21 periods, we would get historical volatility, that is, the actually volatility that occurred over the last 21 trading periods.

We can see the relationship is much closer. The VIX still generally reflects a premium over historical volatility, but big moves in the underlying SPX correlate to big moves in VIX/IV.

Relations


The correlations between these three measures, IV, RV and HV are worth taking a look at. This is using data from January 2000 up to about the end of January 2013:



There is a 91% correlation between historical volatility and VIX. This is higher than the correlation between HV and RV, and HV is really just RV shifted back 21 periods. It seems the relation between VIX and HV is something quite strong.

Now because these are all measures of very similar things, we might expect high correlation. Another metric we could use too look at the relationship is R-Squared.



We can see a relatively strong relationship between historical vol and IV/VIX, at 0.60. The relationship between IV and realized vol is looking a lot weaker.

Conclusion


From all of this, I concluded that the VIX is not a particularly useful indicator of 30 day future realized volatility. However, it is a reasonable estimate of what historical volatility was.

As it turns out, this is not particularly useful, as we can easily calculate historical volatility. It does mean that we can flip things around, and use historical volatility as a reasonable proxy for VIX or market expectations of future realized volatility.

This will turn out to be quite useful, as we may find ourselves trading instruments that do not have nice, readily available VIX equivalents, but that is a topic for a subsequent post.

So ends part 1. These posts turned out very long so I decided to split them up, part 2 should be available here.

The R code for this and the next post is up here.

Thursday, February 14, 2013

OpenGL & openFrameworks Lighting example

I have been wanting to do more with lighting, but did not really feel like I had a reasonable grasp on How It Worked.

I made this app which has the basic lighting types:

  • Ambient Lighting
  • Directional Lighting
  • Point Lighting
  • Spot Lighting

It lets you set the various colours, like the diffuse colour and the specular highlights, and for the material used to colour the sphere you can also set the emissive light.

It is basic but for the basics it is fairly comprehensive, if you aren't quite sure what a diffuse or specular colour is, it will probably be useful. The only thing it doesn't have controls for is light attenuation and I commented out the alpha channel controls as they take up UI space and alpha blending is an app in itself.

There are a couple of oF lighting demo apps but they don't let you easily control the parameters. Mine uses the ofxUI addon so has external dependencies, but should make it a lot easier to see what happens when you shine a blue diffuse light with a red specular highlight on a yellow object.

A screen shot is below. I've made a couple of changes, the biggest is if you are holding down the space bar as you change a color it will change the RGB components together, so you can adjust in terms of black to white.



Source is up here: OpenGL openFrameworks Lighting Demo

These are the instructions from the top of the source:


Spot light defaults to a red light with a blue specular highlight
Point light is similar to a spot light, defaults to green light with white spec highlight
Directional light defaults to blue light with red specular light.
Ambient light defaults to weak grey.
Hold down space when changing a color to adjust black <-> white
Use the 'X Source' toggle to draw the light source
You can also adjust the material properties, shine and the various colours.
Toggles at the bottom of the control panel let you turn lights on/off, and also face culling
It uses ofEasyCam, left click + drag to move around right click + drag to zoom
This app depends on ofxUI and in turn ofXmlSettings

Tuesday, January 29, 2013

Tracking down errors in R


It's that moment we all know and love, somewhere in our code something has gone wrong. We think we have done everything right, but instead of expected glory we find only terse red text lain below our lintel. 

This can be very frustrating, and trouble shooting these issues can often be very time consuming. 

All is not lost. There are a few bits of R that can greatly help finding out what exactly has gone wrong and where, which in turn should suggest a reasonable course of action.

First we will look at some simple methods we can use to track down issues, namely the warn option and traceback, and then we will look at stepping through functions with debug. 

First Steps


I'm going to use an example that has stuck with me from when I first started using R, using neural nets for classification with the iris data.

Let's take a look at the error:

library(nnet)
X <- iris[,1:4]
Y <- iris[,5]
mod <- nnet(X, Y, size=2)

# weights: 19
Error in nnet.default(X, Y, size = 2) : 
  NA/NaN/Inf in foreign function call (arg 2)
In addition: Warning message:
In nnet.default(X, Y, size = 2) : NAs introduced by coercion

Urgh. Really? What the hell does that mean?

We can use the built in traceback() function to see where this error occured

> traceback()
2: nnet.default(X, Y, size = 2)
1: nnet(X, Y, size = 2)

We can see our call to nnet(), which in turn has called nnet.default() and which is where our error has come from. 

In the error output, we can see there was also a warning "NAs introduced by coercion". As we weren't expecting any warnings, let's track down that, as errors tend to compound. 

Warnings


To find out where that message was coming from, we will use options(warn = 2) which will turn warning messages into errors. We can do this by setting the warn option to a specific level, in this case 2. 

The default is warn = 0, which means warnings will be stored until the top level function returns. We could use warn = 1 which wil print the warning as it is encountered, but in this case we want to stop straight away, so will set it to 2.

options(warn = 2)

Let's try again:

> mod <- nnet(X, Y, size=2)
# weights: 19
Error in nnet.default(X, Y, size = 2) : 
  (converted from warning) NAs introduced by coercion
>

Hmm, still coming from nnet.default, let's see if traceback() is offering any new information.

> traceback()
6: doWithOneRestart(return(expr), restart)
5: withOneRestart(expr, restarts[[1L]])
4: withRestarts({
       .Internal(.signalCondition(simpleWarning(msg, call), msg, 
           call))
       .Internal(.dfltWarn(msg, call))
   }, muffleWarning = function() NULL)
3: .signalSimpleWarning("NAs introduced by coercion", quote(nnet.default(X, 
       Y, size = 2)))
2: nnet.default(X, Y, size = 2)
1: nnet(X, Y, size = 2)

We see a whole bunch of extra stuff in the traceback(), but from 3 onwards it appears to be sideffects of having set warn = 2. We do however see our warning has come from nnet.default again, so we will dig into that to see if we can find out what is going on.

Getting dirty with debug()


To do this, we can use the debug function. We will turn on debugging for nnet, which will let us step through the code line by line as it is executed.

> debug(nnet)
> mod <- nnet(X, Y, size=2)
debugging in: nnet(X, Y, size = 2)
debug: UseMethod("nnet")
Browse[2]> 

The Browse> prompt tells us we are in the debugger. The debug: UseMethod("nnet") tells us the next line of code to be executed is UseMethod("nnet"). We could enter 'n' here to continue to the next line, however a convenient default is just hitting enter (i.e. an empty line). 

Browse[2]> 
debugging in: nnet.default(X, Y, size = 2)
debug: {
    net <- NULL
...
    class(net) <- "nnet"
    net
}
Browse[3]> 

Here R has printed out the R source for the function we just entered. We can see we are still at our Browse> prompt, so will continue on by hitting enter again and again:

Browse[3]> 
debug: net <- NULL
Browse[3]> 
debug: x <- as.matrix(x)
Browse[3]> 
debug: y <- as.matrix(y)
Browse[3]> 
debug: if (any(is.na(x))) stop("missing values in 'x'")
Browse[3]> 
debug: NULL
Browse[3]> 

Well this is neat. Hitting enter at the prompt, R shows us each line that is about to be executed. We will continue on hitting enter until we see our error message.

Browse[3]> 
… 
Browse[3]>
debug: if (length(weights) != ntr || any(weights < 0)) stop("invalid weights vector")
Browse[3]> 
debug: NULL
Browse[3]>

One thing of interest, we can see a conditional if statement is about to be run. When the conditional is evaluated as false, meaning the conditioning code won't be executed, we will see the debug: NULL printed out. 

Browse[3]> 
debug: Z <- as.double(cbind(x, y))
Browse[3]> 
Error in nnet.default(X, Y, size = 2) : 
  (converted from warning) NAs introduced by coercion

Well, there is our warning message. Unfortunately we have lost the Browse> prompt, meaning we are no longer inside the function being debugged, but back at the main prompt. 

This is a side effect of our rather aggressive options("warn") setting. Let's tone it down a bit and set it to 1 so the warnings will be printed as they occur, then jump back into debugging.

> options(warn=1)
> mod <- nnet(X, Y, size=2)
debugging in: nnet(X, Y, size = 2)
debug: UseMethod("nnet")
Browse[2]> 
debugging in: nnet.default(X, Y, size = 2)
debug: {
    net <- NULL
Browse[3]> 
Browse[3]> 
debug: Z <- as.double(cbind(x, y))
Browse[3]> 
Warning in nnet.default(X, Y, size = 2) : NAs introduced by coercion
debug: storage.mode(weights) <- "double"
Browse[3]> 

After some time, we get our warning message, and we are still in our debugged function. The warning is coming from the line

Z <- as.double(cbind(x, y))

What could be the problem here? Something is going wrong when nnet is converting x and y to doubles. Let's take a look at them to see if there is anything going on.

We can do this by using get("variable"), where variable is the quoted name of the variable. First let's take a look at other two variables in the function to see how it works:

Browse[3]> get("nout")
[1] 1
Browse[3]> get("ntr")
[1] 150
Browse[3]> 

Looking back at the debug output of each line, we can see they were set as some of the dimension info for X and Y, nout has a value of 1 and ntr has a value of 150. 

Let's take a look at x and y now:

Browse[3]> get("x")
       Sepal.Length Sepal.Width Petal.Length Petal.Width
  [1,] 5.1 3.5 1.4 0.2
  [2,] 4.9 3.0 1.4 0.2
  [3,] 4.7 3.2 1.3 0.2
  [4,] 4.6 3.1 1.5 0.2
Browse[3]>

This is our input X data.

Browse[3]> get("y")
       [,1]        
  [1,] "setosa"    
  [2,] "setosa"    
  [3,] "setosa"    
  [4,] "setosa"    
...
Browse[3]>

And here is our Y data, just as we passed in.

Recall the warning was triggered by this line:

Z <- as.double(cbind(x, y))

Which is converting x and y to the numeric double type. 

We can see that x is numeric data, while y is character data. What happens if you convert character data to numeric? It doesn't seem to make sense, but let's try:

Browse[3]> as.numeric(get("y"))
Warning: NAs introduced by coercion
  [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
 [37] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
 [73] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[109] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[145] NA NA NA NA NA NA
Browse[3]> 

A-ha!

This is just the warning we see, and thinking about it, we can understand why converting strings to numeric data is probably not going to be particularly meaningful.

Let's get out of the debugger and think about what might be going on. Hit 'c' to continue execution 

Browse[3]> c
Error in nnet.default(X, Y, size = 2) : 
  NA/NaN/Inf in foreign function call (arg 2)
>

We are now back at the main prompt. 

Now what??


What is going on? It seems reasonable that we should be able to pass a factor for classification, in fact we are pretty sure that's what we saw being used in the examples in the nnet package documentation.

A careful reading of help(nnet) reveals some details. In particular:

If the response in formula is a factor, an appropriate classification network is constructed ... If the response is not a factor, it is passed on unchanged to nnet.default.

It is possible to pass a factor in Y, but we must use the formula syntax. Looking at the examples, we see the matrix syntax in use as well, however it is transforming the y values using class.ind(). 

Right about now I contemplate calling the police, as I'm pretty sure someone has snuck in and changed the docs while I wasn't looking. There is no way I would make such a simple mistake ...  

Anyway, let's turn off debugging, and see if we can get this working:

> undebug(nnet)
> mod <- nnet(X, class.ind(Y), size=2)
> mod <- nnet(Y~., data=cbind(X, Y), size=2)

Both methods run without error (or warning). Success at last.

Outro


We've seen a few ways we can dig into R and track down where things are going wrong. First is using options(warn=2) to make the R convert warnings to errors, and using traceback() to find out in which function the issue is arising. 

Often, this may be enough to get things back on track, especially if the function causing trouble is small. For more complex issues, we can use debug() which will let us step through the function line by line, inspecting variables and internal state as needed.

It should be said this will only be of use when the source is available to R, i.e. the function in question has been implemented in R and not compiled C/C++/Fortran or whatever that is imported as a shared lib. Deubgging that is more involved and there is some reasonable documentation on this available here.

That's all for now, happy hunting and thanks for stopping by.

Sunday, January 20, 2013

Interactive deJong Attractor

I have been working on this the last few weeks. It is an interactive deJong Attractor, developed using GLSL Shaders and C++ with openFrameworks.

The attractor is generated, then passed as a texture to a shader that builds a curl, then a bunch of particles float around that curl, changing color and opacity based on the attractor output.

It is also connected to a kinect, so the changes that you see on screen are coming from movement captured from the kinect.

I set it up on the dancefloor at a friends party, and included some iPhone footage of it in action at the end of the clip below.

Much of the detail gets lost when looking at the small vimeo embed, so please check out the larger HD version on vimeo here.


Attractor from dizzy pete

I think I have finally gotten deJong attractors and particle fields out of my system, on to the next thing.

Monday, December 17, 2012

Density Plot with ggplot

This is a follow on from the post Using apply sapply and lappy in R.

The dataset we are using was created like so:


m <- matrix(data=cbind(rnorm(30, 0), rnorm(30, 2), rnorm(30, 5)), nrow=30, ncol=3)

Three columns of 30 observations, normally distributed with means of 0, 2 and 5. We want a density plot to compare the distributions of the three columns using ggplot.

First let's give our matrix some column names:

colnames(m) <- c('method1', 'method2', 'method3')
head(m)
#         method1    method2  method3
#[1,]  0.06288358  2.7413567 4.420209
#[2,] -0.11240501  3.4126550 4.827725
#[3,]  0.02467713  1.0868087 4.044101

ggplot has a nice function to display just what we were after geom_density and it's counterpart stat_density which has more examples. 

ggplot likes to work on data frames and we have a matrix, so let's fix that first

df <- as.data.frame(m)
df
#       method1    method2  method3
#1   0.06288358  2.7413567 4.420209
#2  -0.11240501  3.4126550 4.827725
#3   0.02467713  1.0868087 4.044101
#4  -0.73854932 -0.4618973 3.668004

Enter stack


What we would really like is to have our data in 2 columns, where the first column contains the data values, and the second column contains the method name. 

Enter the base function stack, which is a great little function giving just what we need:

dfs <- stack(df)
dfs
#        values     ind
#1   0.06288358 method1
#2  -0.11240501 method1
#…
#88  5.55704736 method3
#89  6.40128267 method3
#90  3.18269138 method3

We can see the values are in one column named values, and the method names (the previous column names) are in the second column named ind. We can confirm they have been turned into a factor as well:

is.factor(dfs[,2])
#[1] TRUE

stack has a partner in crime, unstack, which does the opposite:

unstack(dfs)
#       method1    method2  method3
#1   0.06288358  2.7413567 4.420209
#2  -0.11240501  3.4126550 4.827725
#3   0.02467713  1.0868087 4.044101
#4  -0.73854932 -0.4618973 3.668004

Back to ggplot


So, lets try plot our densities with ggplot:

ggplot(dfs, aes(x=values)) + geom_density()

The first argument is our stacked data frame, and the second is a call to the aes function which tells ggplot the 'values' column should be used on the x-axis.

However, our plot is not quite looking how we wish:


Hmm. 

We want to group the values by each method used. To do this we will use the 'ind' column, and we tell ggplot about this by using aes in the geom_density call:

ggplot(dfs, aes(x=values)) + geom_density(aes(group=ind))


This is getting closer, but it's not easy to tell each one apart. Let's try colour the different methods, based on the ind column in our data frame.

ggplot(dfs, aes(x=values)) + geom_density(aes(group=ind, colour=ind))



Looking better. I'd like to have the density regions stand out some more, so will use fill and an alpha value of 0.3 to make them transparent.

ggplot(dfs, aes(x=values)) + geom_density(aes(group=ind, colour=ind, fill=ind), alpha=0.3)



That is much more in line with what I wanted to see. Note that the alpha argument is passed to geom_density() rather than aes().

That's all for now.




Using apply, sapply, lapply in R


This is an introductory post about using apply, sapply and lapply, best suited for people relatively new to R or unfamiliar with these functions. There is a part 2 coming that will look at density plots with ggplot, but first I thought I would go on a tangent to give some examples of the apply family, as they come up a lot working with R.

I have been comparing three methods on a data set. A sample from the data set was generated, and three different methods were applied to that subset. I wanted to see how their results differed from one another.

I would run my test harness which returned a matrix. The columns values were the metric used for evaluation of each method, and the rows were the results for a given subset. We have three columns, one for each method, and lets say 30 rows, representing 30 different subsets that the three methods were applied to. 

It looked a bit like this

        method1  method2    method3 
[1,] 0.05517714 0.014054038 0.017260447
[2,] 0.08367678 0.003570883 0.004289079
[3,] 0.05274706 0.028629661 0.071323030
[4,] 0.06769936 0.048446559 0.057432519
[5,] 0.06875188 0.019782518 0.080564474 
[6,] 0.04913779 0.100062929 0.102208706

We can simulate this data using rnorm, to create three sets of observations. The first has mean 0, second mean of 2, third of mean of 5, and with 30 rows.

m <- matrix(data=cbind(rnorm(30, 0), rnorm(30, 2), rnorm(30, 5)), nrow=30, ncol=3)

Apply


When do we use apply? When we have some structured blob of data that we wish to perform operations on. Here structured means in some form of matrix. The operations may be informational, or perhaps transforming, subsetting, whatever to the data.

As a commenter pointed out, if you are using a data frame the data types must all be the same otherwise they will be subjected to type conversion. This may or may not be what you want, if the data frame has string/character data as well as numeric data, the numeric data will be converted to strings/characters and numerical operations will probably not give what you expected. 

Needless to say such circumstances arise quite frequently when working in R, so spending some time getting familiar with apply can be a great boon to our productivity.

Which actual apply function and which specific incantion is required depends on your data, the function you wish to use, and what you want the end result to look like. Hopefully the right choice should be a bit clearer by the end of these examples.

First I want to make sure I created that matrix correctly, three columns each with a mean 0, 2 and 5 respectively. We can use apply and the base mean function to check this. 

We tell apply to traverse row wise or column wise by the second argument. In this case we expect to get three numbers at the end, the mean value for each column, so tell apply to work along columns by passing 2 as the second argument. But let's do it wrong for the point of illustration:

apply(m, 1, mean)
# [1] 2.408150 2.709325 1.718529 0.822519 2.693614 2.259044 1.849530 2.544685 2.957950 2.219874
#[11] 2.582011 2.471938 2.015625 2.101832 2.189781 2.319142 2.504821 2.203066 2.280550 2.401297
#[21] 2.312254 1.833903 1.900122 2.427002 2.426869 1.890895 2.515842 2.363085 3.049760 2.027570

Passing a 1 in the second argument, we get 30 values back, giving the mean of each row. Not the three numbers we were expecting, try again. 

apply(m, 2, mean)
#[1] -0.02664418  1.95812458  4.86857792

Great. We can see the mean of each column is roughly 0, 2, and 5 as we expected.

Our own functions


Let's say I see that negative number and realise I wanted to only look at positive values. Let's see how many negative numbers each column has, using apply again:

apply(m, 2, function(x) length(x[x<0]))
#[1] 14  1  0

So 14 negative values in column one, 1 negative value in column two, and none in column three. More or less what we would expect for three normal distributions with the given means and sd of 1. 

Here we have used a simple function we defined in the call to apply, rather than some built in function. Note we did not specify a return value for our function. R will magically return the last evaluated value. The actual function is using subsetting to extract all the elements in x that are less than 0, and then counting how many are left are using length

The function takes one argument, which I have arbitrarily called x. In this case x will be a single column of the matrix. Is it a 1 column matrix or a just a vector? Let's have a look:

apply(m, 2, function(x) is.matrix(x))
#[1] FALSE FALSE FALSE

Not a matrix. Here the function definition is not required, we could instead just pass the is.matrix function, as it only takes one argument and has already been wrapped up in a function for us. Let's check they are vectors as we might expect.

apply(m, 2, is.vector)
#[1] TRUE TRUE TRUE

Why then did we need to wrap up our length function? When we want to define our own handling function for apply, we must at a minimum give a name to the incoming data, so we can use it in our function.

apply(m, 2, length(x[x<0]))
#Error in match.fun(FUN) : object 'x' not found

We are referring to some value x in the function, but R does not know where that is and so gives us an error. There are other forces at play here, but for simplicity just remember to wrap any code up in a function. For example, let's look at the mean value of only the positive values:

apply(m, 2, function(x) mean(x[x>0]))
#[1] 0.4466368 2.0415736 4.8685779

Using sapply and lapply


These two functions work in a similar way, traversing over a set of data like a list or vector, and calling the specified function for each item.

Sometimes we require traversal of our data in a less than linear way. Say we wanted to compare the current observation with the value 5 periods before it. Use can probably use rollapply for this (via quantmod), but a quick and dirty way is to run sapply or lapply passing a set of index values.

Here we will use sapply, which works on a list or vector of data. 

sapply(1:3, function(x) x^2)
#[1] 1 4 9

lapply is very similar, however it will return a list rather than a vector:

lapply(1:3, function(x) x^2)
#[[1]]
#[1] 1
#
#[[2]]
#[1] 4
#
#[[3]]
#[1] 9

Passing simplify=FALSE to sapply will also give you a list:

sapply(1:3, function(x) x^2, simplify=F)
#[[1]]
#[1] 1
#
#[[2]]
#[1] 4
#
#[[3]]
#[1] 9

And you can use unlist with lapply to get a vector.

unlist(lapply(1:3, function(x) x^2))
#[1] 1 4 9

However the behviour is not as clean when things have names, so best to use sapply or lapply as makes sense for your data and what you want to receive back. If you want a list returned, use lapply. If you want a vector, use sapply.

Dirty Deeds


Anyway, a cheap trick is to pass sapply a vector of indexes and write your function making some assumptions about the structure of the underlying data. Let's look at our mean example again:

sapply(1:3, function(x) mean(m[,x]))
[1] -0.02664418  1.95812458  4.86857792

We pass the column indexes (1,2,3) to our function, which assumes some variable m has our data. Fine for quickies but not very nice, and will likely turn into a maintainability bomb down the line. 

We can neaten things up a bit by passing our data in an argument to our function, and using the special argument which all the apply functions have for passing extra arguments:

sapply(1:3, function(x, y) mean(y[,x]), y=m)
#[1] -0.02664418  1.95812458  4.86857792

This time, our function has 2 arguments, x and y. The x variable will be as it was before, whatever sapply is currently going through. The y variable we will pass using the optional arguments to sapply

In this case we have passed in m, explicitly naming the y argument in the sapply call. Not strictly necessary but it makes for easier to read & maintain code. The y value will be the same for each call sapply makes to our function. 

I don't really recommend passing the index arguments like this, it is error prone and can be quite confusing to others reading your code. 

I hope you found these examples helpful. Please check out part 2 where we create a density plot of the values in our matrix.

If you are working with R, I have found this book very useful day-to-day R Cookbook (O'Reilly Cookbooks)