Coroutintes for simulations

15 May 2010 18:00

Returning to the problem of simulating investment strategies, I was struggling with Clojure. I have spent a lot more time with it and even written a basic asynchronous web-server with it, but there are still things which don't feel natural to me.

The task of simulating an investment strategy seemed ill-suited. After each simulated day, my simulation needs to update its ideas about the world and decide whether to invest based on that. Whilst it could re-compute everything for each day from all the previous history, keeping state and updating that state based only on the current day seemed a more sensible option.

With Clojure, the above would probably require something which could pass some functions their state and the next day's data, and then have the functions return their new state and a function deciding whether to invest. (This is a function of the following day's opening price).

What feels more natural is to have the state stay in place and to just pass in the next day's data and have the investment decision returned. This sounded like a perfect job for coroutines, so I decided that it was time to switch back to Python. There is an excellent tutorial on using coroutines in Python (pdf).

Here are some of the coroutines which I ended up creating:

import functools

def coroutine(f):
@functools.wraps(f)
def start(*args,**kwargs):
cr = f(*args,**kwargs)
cr.next()
return cr
return start

def mean(l):
if l:
return sum(l)/len(l)

@coroutine
def simpleMovingAverage(days):
cached = []
while 1:
cached.append((yield mean(cached)))
if len(cached) > days:
cached.pop(0)

@coroutine
def exponentialMovingAverage(f):
a = None
while 1:
if a is None:
a = yield
a = f*(yield a) + (1-f)*a

@coroutine
def difference(c1, c2):
v = yield
while 1:
v = yield c1.send(v) - c2.send(v)

@coroutine
def isPositive(c):
v = yield
while 1:
v = yield c.send(v) > 0

Combining these coroutines enables complex things to be expressed with good encapsulation. For example, suppose I want a coroutine which returns True iff the 10 period simple moving average is above the 30 period exponential moving average. Then I can combine them like this:

trend = isPositive(difference(simpleMovingAverage(10),
exponentialMovingAverage(2.0/31)))

With this it is easy to simulate what would happen if I held stocks when the short-term moving average is above the long-term moving average.