Fitting a mixture of independent Poisson distributions

This is an example from Zucchini & MacDonald’s book on Hidden Markov Models for Time Series (exercise 1.3).

The data is annual counts of earthquakes of magnitude 7 or greater, which exhibits both overdispersion for a Poisson (where the mean should equal the variance) as well as serial dependence.

The aim is to fit a mixture of m independent Poisson distributions to this data, using the non-linear minimizer nlm to minimize the negative log-likelihood of the Poisson mixture.

Sounds easy right? I have not done much optimisation stuff so it took me longer than it probably should have. There’s little better for the ego then getting stuck on things by page 10.

Constrained to unconstrained 

There are two sets of parameters, the lambda parameters to the Poisson distributions and the delta mixing distribution, the latter giving the proportions of each distribution to mix.

The Poisson parameters must be greater than zero and there are m of them. The mixing distribution values must all sum to one, and the first mixing value is implied by the m-1 subsequent values.

These are the constraints (lambdas greater than zero, deltas sum to one), and there are 2m – 1 values (m lambdas, m-1 deltas). We need to transform these to unconstrained values for use with nlm (becoming eta and tau respectively).

The formulas for transformation are reproduced from the book here

These transformed values are combined in to one parameter vector and passed to nlm. We also pass a function that calculates the negative log likelihood.

In the function, we first convert the parameters from the unconstrained “working” form to the constrained “natural” form using the inverse transform, and then use these values to do the likelihood calculation.

You can see the code here.

Outro

Now that I look over the code and write this out it does seem all fairly straightforward and relatively trivial, but it was new to me and took a while, so I figured other people might find the example useful as well.

The code has three examples for the cases m = 2, 3, and 4. Obviously feel free to mess around with the initial parameter estimates to see what comes up. 

It does seem to match the fitted models from the book, which I have copied here:

For all the details, check out the book, chapter 1!

A quick look at FX realized vol

Much has been said about the decline in volatility. At the moment I am very active in FX spot trading and as a generalization do better the more vol there is.

 I wanted to see how things stood on the crosses I am most active in, namely EUR/USD, GBP/USD and USD/JPY.

 I took hourly data from FxPro (not my broker, nor an endorsement), calculated volatility as the high minus the low, and summed the total for each day. You can think of it as how many pips were on offer if one could correctly call the high and low of each hour of each day.

 All up there is about 90 days of data, so it covers roughly the last four months. I also took the average of the last five days which are the red X’s on the box plots. We are here.

As you can see, vol is below average. It was quiet week overall, bank holidays in the US and parts of Europe, and only a moderate amount of data coming out. Next week should be a bit busier I think.

Since I had all the data I also took at look at the average hourly RV per day.

I don’t want to read too much into this chart but things have been quiet. I read somewhere else fx vol is approaching levels of 2007 which was a very quiet time indeed.

Some R code is up here, data is here.

Fractional sums of Perlin Noise

I wanted to mess around with fractional sums of Perlin noise, so made a little openFrameworks app to better understand what is going on.

Frequency

Frequency comes up quite a bit in the following discussion, and more or less means the rate at which something goes up and down. All you really need to keep in mind is:

High frequency noise goes up and down quickly. It looks like this (where black is 0 and white is 1):

Low frequency noise goes up and down slowly.  It looks like this:

You can see low frequency noise takes a while to get from zero to one and vice versa, while high frequency noise does it much more often.

ofNoise Inputs and Outputs

In general it’s a good starting point to pass normalized coordinates to the noise function.

You can use ofNoise() to get values ranging from [0, 1], or ofSignedNoise() to get values from [-1, 1].

There are a bunch of functions for various dimensions as well.

Fractional Sums

Low frequency noise will give a nice undulating look, but often it is boring. High frequency noise is more interesting, but can be a bit too chaotic. What we want is a nice combination of both.

Using an example from Paul Bourke, we take several instances of noise at increasing frequency and combine them to get the effect we are after.

To do this there are three parameters, an octave count, alpha, and beta.

The octave count is how many layers of noise we will be adding together. This will typically range from 1 to 8.

Each layer of noise is generated at a higher frequency than the one before, which is where the name octave comes from. 

Beta controls the frequency of noise, the larger it is, the higher the frequency of the noise, i.e. the faster it goes up and down.

Adding these together works, but you might find the higher frequency noise overpowers the lower frequency noise.

This is where the third parameter, alpha, comes in, controlling how much of the noise from the current octave ends up in the final sum.

In rough code it looks like this:

double sum = 0.0;
for (int n  = 0; n < octaveCount; n++) 
            sum += 1/alpha^n * noise(n * beta * x, n * beta * y);

Lets say alpha and beta are both 2. In the call to noise(), the n * beta term will get progressively larger for each successive layer of noise. We are increasing the frequency of the generated noise.

However, with alpha as 2, we are adding successively less of each octave, as we get 1/2, 1/4, 1/8, 1/16 … reducing the magnitude of the higher frequency terms.

Taken together, we see the higher frequency noise contributes less and less to the final sum. We end up with a nice smooth but varying noise map to use for displacement or whatever we want.

Normalizing

Summing several octaves means you will typically get values greater than 1 (and less than -1 if you are using signed noise).

There’s a few ways you can normalize these, take a look at the commented out code. You might not want to map the minimum to zero as it can cause jumps when the lower bound of the summed noise changes.

The app 

An app to play with is up here

The small images show the individual octaves, and the big image is the final result.

Also I used the noise function from that app to do vertex displacement of a mesh and put a clip up here:

That’s it for now. Later on I will give an example of making seamless noise loops. I have a bunch of other stuff going on which I will write about when they are finished off, but I am really looking forward to sharing.

As always you can find me being rude and unprofessional on twitter. I love hearing from people, so tweet me pics of your rad noise stuff. 

Cheap Tricks - Interactive Dance Floor Application with Kinect and openFrameworks

I took a short break from a longer term project I am working on to make something for a New Years Eve party dance floor. I didn’t have a lot of time and it shows one of my favorite aspects of working with openFrameworks: you get a lot of stuff for free via the addons available, making it very quick to implement ideas.

You can see a little demo of the project here:


Cheap Tricks from dizzy pete.

It uses a Kinect and OpenCV to get the blobs (in this case people and my Ikea furniture) then turns them into triangles that are audio reactive, as well as adding in some visual effects with no real purpose other than to look cool.

Kinect Input

Reading from the Kinect is relatively straightforward, and more or less lifted straight from Kinect examples that come with oF. It uses the depth information in a grayscale OpenCV image and does some simple thresholding based on depth settings.

This grayscale image is passed to the stock OpenCV blob finder, then the blob information is passed to ofxTriangle to get a set of triangles representing the blob. I initially was using ofxDelaunay but found ofxTriangle just worked a little better for what I was after.

Audio

The audio reactive part is done via ofxFFT by julapy (there are a few out there so make sure you have the right one if you’re taking a look). I originally wanted to do beat detection and have the figures pulse in time with the beat, but that was a little too involved given the time constraint, so decided just to stick with the very low end.

I experimented with various numbers of FFT bins but ended up using the default buffer size of 512, bin size 256 and averaging the lowest three bins. This is the variable ‘low’ that is used throughout the code.

Ideally I would have had a proper audio device feeding in from the DJ’s mixer, but as it turned out I was using the default microphone on my MacBook for audio input, which worked much better than I expected.

Based on the current amount of low-end audio detected, each vertex of each triangle is displaced from its center by some amount. Triangle sizes vary; the torso of the detected body has larger triangles than the hand or arm, how much displacement to be applied is a function of the triangle area.

To be precise it is the square root of the triangle’s area, and can be tweaked by a parameter called jitterSize in the gui. I don’t actually displace each vertex, just two out of three. The code to do all third is there but commented out, I just liked the way it looked doing only two.

Once the vertex displacement has been calculated, the triangle is added to a mesh along with texture coordinates, which are used for coloring. But more on that in a second.

The basic triangle form

Background

By this stage the basic audio-reactive triangle people thing I was after was done, but it still looked a bit flat.

Poking around some OpenCV examples, I saw examples of motion detection that looked kinda cool.

First it takes the difference between the current frame and last frame, which gives areas where motion has been detected. This is then scaled up as a function of low-end audio and GUI parameter glowExtra, then added to a buffer image. The buffer image gets decayed over time, so things fade out as the motion fades.

This looked ok, but I wanted it to be a bit smoother so used ofxFboBlur to add in excessive amounts of Gaussian blur to the motion buffer image. I had to change its source a bit just to play a bit nicer with FBO sizes and drawing, which is why it appears in the source folder and not as an addon.

The blurred motion detection background

Color

Coloring was probably the hardest part for me, and even then it was pretty easy.

There are lots of different ways this could be done, but I settled on using a two-color gradient rendered to an FBO, using a mesh with vertex colors. Each frame, the mesh vertex colors are updated by cycling through the Hue component in HSB space, then the updated mesh is rendered to the FBO.

This FBO gets used in two places. First it is blended with the grayscale blurred motion detection background using a multiplicative blend mode, which takes care of coloring in the swishy aura type effect in the background.

Then the FBO texture reference is bound and the mesh containing our displaced triangularized blobs is drawn as a wireframe. This is where having texture coordinates for each vertex in the mesh is important, though fiddling around I noticed I don’t need tex coords if ofBackground() has been called. Who knows.

Regardless, when the wireframe is drawn, it is colored based on the bound FBO texture reference, which takes care of coloring in our triangle people.

The two layers colored and drawn together

Outro

So that is just a quick run through of how it works. I used a lot of addons, and all up it took me about 2 days. I started on the evening of the 29th and was finished by the afternoon of the 31st, just in time. Without the oF addons it would’ve taken me well over a week, possibly many weeks.

People at the party really seemed to enjoy it as well, which is always a great feeling.

I used soundflower during development (and oF 0.7.3 eep), but updated it to use oF 0.8, which is why you see ofxSoundStream in the code.

If you don’t have a Kinect it will use the device camera, but it does not work as well. You can use the far threshold setting for the OpenCV threshold in that case.

The source is up here.

Also I have put the source for Particles up here as well

You can find me on twitter here !


Particles 1 from dizzy pete.

A volatility filter using historical vol

We have been looking at a way to improve risk adjusted returns by using a volatility filter. Although we could use VIX or equivalent, it turns out that historical volatility will work just as well, if not a little better.

You can see part 1 here Digging into the VIX, and part 2 here What can we use VIX for?

Although the mean return of how we slice things is zero, the distribution of returns is wider for higher readings of our relative measure of volatility. High volatility begets high volatility, at least for our purposes.

By staying out during periods of higher relative volatility, we aim to reduce drawdowns and the volatility of our returns, leading to better risk adjusted results.

A plus of using HV over some external measure like VIX is that it is readily available for any underlying. This means such a filtering technique can be applied to whatever it is we are trading.

Performance

Below is a table with two comparisons, the first compares the HV filter to buy and hold.  Although performance is generally better, we still get some pretty big drawdowns.

The second adds in a 200 period moving average, which is a reasonably strong way of protecting against downside. Again we can see lower volatility and smaller drawdowns with the addition of a vol filter.

I used a 3 month/63 day look back for our relative volatility measure. I haven’t really dug in to what happens when volatility remains elevated for extended periods of time.

I also did not experiment much with the threshold for where we draw the line on ‘high’ relative volatility. I use 0.6 as the cutoff because I originally split things into quintiles when making the first charts.

I also ran this for RUT and NDX over the same period

These results are all frictionless, don’t factor in dividends, return on cash, etc, etc. I don’t consider this viable as a standalone indicator, but something that can be used along side other factors like rotational strategies, or as a potential tool if you are looking for lower volatility.

The source is up here, feel free to have a play around with it and see how you go.

Thanks for reading, 'till next time.

Another cut at market randomness

I have some background in computer security and one day found myself tasked with assessing the quality of randomness for session id tokens generated by popular web frameworks (namely Java and .NET). As it turns out, NIST have developed a series of tests for just this purpose detailed here.

As a non-believer in the absolute randomness of markets, I thought I might take a look at the series of returns from SPX/GSPC to see how they held up.

About the tests

The NIST tests are designed to be run on bit streams from a random number generator, so some liberties had to be taken. I am mainly interested in the predictability of up days vs down days, so I first generated a series of returns using quantmod then encoded them using a 1 for an up day and -1 for a down day and run the tests on the encoded stream. The NIST tests use -1 in place of a 0.

The main test is called the monobit test, to quote from the document:

The focus of the test is the proportion of zeroes and ones for the entire sequence.  The purpose of this test is to determine whether the number of ones and zeros in a sequence are approximately the same as would be expected for a truly random sequence.  The test assesses the closeness of the fraction of ones to ½, that is, the number of ones and zeroes in a sequence should be about the same.  

The tests are done using a significance level of 0.01, so come with a good degree in confidence assuming the underlying method is sound. 

One caveat is the length of the runs compared and how it relates to the distributions used to model the results. For the monobit test, the suggested input size is 100, requiring 101 days of data to determine 100 up or down days. If we were looking at 32 bit integers, 100 bits would only be 3 "full" random numbers, so arguably we would want to look at shorter time periods (e.g. 3-5 days of data). Given the difficulties around distributions which require a large n, I thought I would vary the significance level instead of a lower n, as our requirements are not as stringent as those for cryptographic random numbers.   

Results

At a basic level, this series does appear to be random, at least the vast majority of the time with n = 100 and alpha = 0.01. My confirmation bias was very upset with this. 

However, if we plot the proportion of runs deemed random vs the significance level, we see the proportion rising as one might expect. One thing that remains unexplained is why this appears to rise in steps rather than something more linear, though I expect this to be a side affect of either methodology or the normalisation done by the tests. I also took a look at the weekly data, which tends to a greater proportion of non-random runs quicker than daily data.

I am interested in the applications of machine learning to financial markets. Close to close returns that we have been looking at here are not the only information we have available, nor are they what I trade on a personal level. Also this is only one price series, and one could argue that in practise it is not actually a tradable series. 

Close to close returns are very useful in lots of applications, but if we are trying to build some predictive model we might need to look for more predictable pastures. Machine learning algorithms are great, but can only do so much. Finding some better potential inputs is what I will take a look at next.

The test code is available on github here: R Monobit test. Would be very interested to hear if anyone else takes a look.

I also took a visual and binomial look at randomness of the series in this post: A visual look at market randomness.

Oh and in case you were wondering, the web session tokens all turned out to be very strong.

Mebane Faber Tactical Asset Allocation in R

In 2006 Mebane Faber published a great piece of research detailing an asset allocation system that was both very easy to understand and implement, as well as carrying very respectable risk adjusted returns.

The details are available in his paper A Quantitative Approach to Tactical Asset Allocation and were further expanded on in his book The Ivy Portfolio both of which are must reads.

The short version is to use diversified asset classes, long only, and only long when the price is above the 10 month simple moving average (approx 200 day). The assets he tests are U.S. Stocks, International Stocks, U.S. Government Bonds, Commodities and Real Estate, accessible via ETFs.

A rotational extension can also be added by investing only in the top 1-3 asset classes showing some degree of relative strength, which is defined as the average of 3, 6 and 12 month returns. They must also be over the 10 month SMA to be candidates.

The system updates monthly at the end of the month, it is about as hands off as you can get for active management.

There is an ETF for those so inclined, GTAA, but I am experimenting with a put selling implementation, which I might start tracking here month to month. I wrote a small R script using quantmod to display the relevant information for given symbols, which should be available here: Tactical Asset Allocation R script

The output looks like this:


  Sym         R3m         R6m        R12m  Close     AvgRet OverMA
4 VNQ  0.09295631  0.22412597  0.08488552  63.65 0.13398927   TRUE
1 VTI  0.11671109  0.22466699  0.05037598  72.26 0.13058469   TRUE
2 VEU  0.10908623  0.13282091 -0.10915250  44.22 0.04425155   TRUE
5 DBC  0.07048208  0.11194076 -0.05767911  28.80 0.04158124   TRUE
3 IEF -0.02193049 -0.01718305  0.10473673 103.28 0.02187440   TRUE


Let me know if you have any comments or find any bugs. 

Machine Learning Examples in R

This is a post that has been a long time in the making. Following on from the excellent Stanford Machine Learning Course I have made examples of the main algorithms covered in R.

We have Linear Regression

Followed by Neural Networks

And Support Vector Machines

One remaining item is Logistic Regression, I am yet to find a library in R that behaves as I want, so that will come at some future date. I've been sitting on this post for ages and got sick of waiting. As an aside I find the documentation in R to be variable at best, which can make it somewhat of a pain to work with. When it is good, yes it can be very good but often it is quite poor ...

R is great for data analysis and exploration, but I have found myself moving back to python for many more mundane tasks.

Anyway for those interested in the code, I have put it on Github. The data is from an exercise in the Stanford course, and by tweaking the parameters I really got a good feel for how the various algorithms work in practise.

Once I finish my backtesting engine I will probably put it up on Github as well, and then I can start digging into the applications of ML techniques for trading systems.