Wednesday, June 15, 2016

Domino Data Lab python/bash hacks

In a prior post I covered the data science platform provided by Domino Data Labs. I'm still a fan and here are a few minor gimmicks I've found useful when working with the product, starting with the trivial. You may prefer to cut and paste from domino_data_hacks on github.

(A run utilizing 8320 processors on Amazon. See below for a bash script to do this.)

Starting with the basics... Well I said I would start with the trivial. For readability:

import os
def domino_run_id():
        return os.environ['DOMINO_RUN_ID']
        return None

def running_on_domino():
    return ( domino_run_id() is not None )

def running_on_local():
    return not running_on_domino()
Domino environment variables While we are at it:
current_project       = os.environ['DOMINO_PROJECT_NAME']
current_project_owner = os.environ['DOMINO_PROJECT_OWNER']
current_run_id        = os.environ['DOMINO_RUN_ID']
current_run_number    = os.environ['DOMINO_RUN_NUMBER']
domino_project_path   = '{0}/{1}'.format(current_project_owner, current_project)

Shelling out

Due to vagaries I don't fully understand but are undoubtedly related to security, permissions on files might not be what you expect. To work around this you might go so far as to set the permission right before you need it:["chmod", "+x", ""])["chmod", "+x", "stockfish_"+binary])
    cmd =  ' '.join( ['./' ,fen, str(seconds) , binary, str(threads), str(memory) ] )
    print cmd cmd, shell=True )
Now unfortunately if your code is not in the same project as your data (which I recommend, see below) this won't work, at least if run from the data project. One can make a temporary copy of the script you are shelling out to and also call that just before shelling out.
mkdir -p "${project}/data_project/bash_mirror"
mkdir -p "${project}/data_project/bash_mirror/interface"
cp -R ${interface_dir}/* ${interface_mirror}
Yeah, not the most elegant but it works. There is probably a better way. One of the domino engineers suggested adding sh in front of commands (see this note) but whatever you try, be aware that in Domino, file permissions set on a run in one project will not be preserved when you import that project into another.

Setting paths in bash

For the times when you don't want to rely on Domino variables explicitly

if [[ $OSTYPE == *"darwin" ]]  # Your machine -ne AWS hopefully :) 
    # We're local
    default_project="/Users/projects/my_project"   # Hardwired
    # We're on AWS

# allow override

# Then do something...       
Consistent path names Incidentally a project importing other projects can use full paths such as
but a project that does not share any other projects cannot. To avoid inconsistency, just make sure every project imports one other project, even if it is a dummy project.

Drop-in multiple endpoint functions

Domino allows only one endpoint function per project. To allow the client to call any functions you care to drop in your endpoint file instead (by their names), include these four lines of code at the top of the same and register dispatcher as the official endpoint function.

def dispatcher( func, *args ):
    """ As Domino allows one endpoint per project """
    module  = sys.modules[ __name__ ]
    endpoint_func = getattr( module, func )
    return endpoint_func( *args )
The price you pay is specifying the function you really want as the first parameter in any call. This is not recommended for national security applications.

Safe sync

Very, very occasionally a failed sync can leave client and server in a state where it is inconvenient or confusing to revert to a previous code version. If you're a nervous Nellie like me there is a simple way to reduce the chance of any code confusion:

  1. Copy your source code to a second domino project
  2. Sync both with two separate domino syncs
Long odds against both sync's failing at all, much less in a confusing way. Obviously this needs to be scripted (and turned into a one-click as below) or you'll never actually do it. I've provided a little script that you might wish to modify to your taste. It does even more, manually saving as many previous versions as you care for. It relies on a slightly hacky script which again, might be swapped out for a "better" way as you see fit, such as integration with your source control. Again, not really an issue unless you commingle code with huge piles of data. One click sync Moving on to profound observations that don't really have anything to do with Domino per se, I'll point out that you'll probably do the backup/sync operation so frequently you'll end up trying to hide a terminal window in the corner of your screen. On a mac we can do a little better, creating a one click app launched from the dock.
  1. Chmod your backup/sync script so it is executable
  2. Rename it with .app extension, so that you can then drag it into the mac dock
  3. Rename back to original .sh extension
  4. In finder, right click on "Get Info"and under "Open with" menu, select the terminal application
Unfortunately you'll have to repeat the last step every time you edit this file as your IDE will likely revert it - but that shouldn't be necessary too often.


To use the server's state use domino reset. To use the local state use domino restore. More on recovery of larger projects below.

Bash script to run job on Domino and wait for completion

The script uses the Domino API to start a run and then wait for it to complete.

# Run a job on Domino and wait for it to finish
# Usage
# /mnt/USER/MYPROJECT/ my_arg1 my_arg2


# Send request to start job to domino
curl -X POST \ \
-H 'X-Domino-Api-Key: 7DunTrump1BcouldCnTZVhitcaKBarnQDooRwithAWrifle' \
-H "Content-Type: application/json" \
-d '{"command": ["'"${cmd}"'", "'"${arg1}"'", "'"${arg2}"'"], "isDirect": false}' > ${temporary_response_file}
echo "Sent command to start jobs with arg1 ${arg1} and arg2 ${arg2}."
runId_quoted=$(grep -oE '"runId":"(.*)",' ${temporary_response_file} | cut -d: -f2)
echo "The runId is ${runId}"
rm ${temporary_response_file}

# Now poll until done
while true;
    sleep 60s
    echo "Polling ..."
    response=$(curl$runId \
    -H 'X-Domino-Api-Key: 7DuVOtdgeforzhimFanDitYouWEFfault' | grep -oE '"isCompleted":.*')
    if [[ ${response} == *true* ]]
        echo "Job $runId has finished"

Aside: passing bash variables in curl requests with "'"${}"'"

Incidentally, am I the only one who found this slightly troublesome? Here's ten minutes of my life I'm donating to you:

curl -X POST \ \
-H 'X-Domino-Api-Key: QAYeqRTrump2is7aLdangerousMVUSClownVMSF' \
-H "Content-Type: application/json" \
-d '{"command": ["", "'"${key}"'"], "isDirect": false}'
Note the Double-Single-Double quoting of the bash variable.

Lazy man's map

Want to use more than one machine at once? I discovered that the easiest way is have each job write to the main branch. Domino sync will take care of the syncing that way, and you need only filter each run by some key (say {A..Z}).

!/usr/bin/env bash

# Start multiple jobs on Amazon

sz=${1:-full}      # Just an example of a parameter passed to all jobs 

# Pre-requisite is a file with one line containing space separated keys which break up the jobs - see the bach hack above

read -a kys <<<$(head -n 1 ${ordering_file})

for ky in ${kys[@]}    # or just {A..Z} if you don't care about ordering
    sleep 10s    # Give AWS a little time so jobs get ordered the way you expect
    curl -X POST \ \
    -H 'X-Domino-Api-Key: 7DunjmdTrumph1BwillTZdeeKstoyaVQtheWVnCountry3m4F' \
    -H "Content-Type: application/json" \
    -d '{"command": ["", "/mnt/YOUR_USER", "'"${sz}"'", "'"${ky}"'"], "isDirect": false}'
    echo "Set command to start jobs with source_filter ${ky}"

Ordering your jobs so large input files go first

The bash script provides an easy way to divide up your input data by letter {A..Z} say and order them by size. Hack as you see fit.

#!/usr/bin/env bash

# Create a file which sorts data sizes by letter, so we can order the // jobs sensibly


# Delete the old statistics file if it exists
if [[ -e ${config_dir}/${tmp_file} ]]
   rm ${config_dir}/${tmp_file}

# Create a file with one line per letter:
#      Size    key
#      1223411 A
#      1231233 B
for x in {A..Z}
   data_sz=$(du -c ${data_dir}/${x}* | awk '/./{line=$0} END{print $1}')
   echo "${data_sz} ${x}" >> ${config_dir}/${tmp_file}

# Sort, extract the letters, and convert to a single row
sort -r -t " " -k 1 -g ${config_dir}/${tmp_file} | awk '{print $2}' | tr '\n' ' ' > ${config_dir}/${ordering_file}

rm ${config_dir}/${tmp_file}

# Now you can easily read the ordered keys into a bash array:
#        read -a ordering <<<$(head -n 1 ${config_dir}/${ordering_file})
Adjust as you see fit. I use this in conjunction with the previous hack to reduce wall clock time of big jobs.

Lazy man's map-reduce

You'll often want to wait for the results to come back so you can get on with dependent tasks. I've posted a polling version of the task launcher as well. To poll directly from bash you can do something like this:

while true;
    sleep 1m
    for runId in ${runIds[@]}
        if [[ ${finished_runs} == *"$runId"* ]]
            : # No need to check again after it has finished
            response=$(curl$runId \
            -H 'X-Domino-Api-Key: 7DunjmIbhopeWKTrumVpgDieszAFterrWXibleWQEdeath' | grep -oE '"isCompleted":.*')
            if [[ ${response} == *false* ]]
            elif [[ ${response} == *true* ]]
                echo "Job $runId has finished"
                finished_runs="$finished_runs $runId"
                echo "We have a problem - did not expect this"

    if [[ ${status} == "finished" ]]
        echo "All finished"
# Now do the reduce step because you have all the data ready
Thanks Jonathan Schaller for helping me fix this.

Pull runId out of JSON response

As a minor digression ... there are better general purpose JSON parsers like jq, but this is good enough. Send the curl response to ${response_file} and then:

    runId_quoted=$(grep -oE '"runId":"(.*)",' ${response_file} | cut -d: -f2)
Grant one project write access to another Project A imports project B's files. So A can easily read B's files. However, A cannot write to B. The way I hack around this is to share A's files with B, and then in project B create a script to do the copying of A's files over.
echo "Chilling for a minute so Project A can finish syncing"
sleep 1m
cp -r /mnt/${USER}/Project_A/data/results_${key}* /mnt/${USER}/Project_B/data && echo "Success"
ls -l /mnt/${USER}/Project_B/data/results_${key}*
Let's suppose the above sits in Project B with the name In order to drive this from Project A I write a little "beam me up Scotty" script:
#!/usr/bin/env bash
curl -X POST \$USER/mapreduce/runs \
-H 'X-Domino-Api-Key: QqTYrumPoisAlmp2NutterisWEVtruly' \
-H "Content-Type: application/json" \
-d '{"command": ["", "'"${key}"'"], "isDirect": false}'
This uses the Domino API to kick start the process. This is a kludge and only really works at the end of a run because Project A's files (assuming they have changed) can't be seen by Project B until after the sync occurs. A more sophisticated approach (which also solves other issues) is to wrap the domino task, for example in a luigi task as in this example by, you guessed, Jonathan Schaller. As noted above, you don't necessarily need an explicit collection step if you write results to the main branch.

Domino API - job status

Checking on a job...

domino_api            = Domino(domino_project_path)

def get_run_info(run_id):
    for run_info in domino_api.runs_list()['data']:
        if run_info['id'] == run_id:
            return run_info
Thanks Jonathan.

Domino API - data status

def get_blob_key(commit_id, filepath):
    dir, filename = os.path.split(filepath)
    files = domino_api.files_list(commit_id, path=dir)['data']
    for file in files:
        if file['path'] == filepath:
            return file['key']
def file_exists_in_commit(filepath, commit_id):
    # filepath is relative to project root
    files_list = domino_api.files_list(commit_id, path='/')['data']
    for file in files_list:
        if filepath == file['path']: return True
    return False
Thanks again.

Rich man's map-reduce

As noted above, another way to handle large pipelines is to mixin a Domino task into some pipeline framework, such as Luigi. See Jonathan's example at here

Ignore it and it will go away

I assume the reader is familiar with the usage of .dominoignore. But a classic catch-22 arises when you forget to ignore files. Too many are generated and syncing is a pain - including syncing of the new .dominoignore file that would get you out of the predicament. Sometimes I have had difficulty editing the remote version of .dominoignore through the web browser - which is the obvious solution - and that makes it really difficult to wiggle out of the space issue. So I just create a launcher for it.

#!/usr/bin/env bash
echo $1 >> /mnt/USER/PROJECT/.dominoignore
Trivial but exceedingly useful.

Maintaining a high level of motivation

Just had a nine hour experiment crash at the final stage? Remember that we all face challenges, and none more difficult than those solved by

Wednesday, April 20, 2016

Filtering bond and credit default swap markets

I presented this talk at the NYU-Tandon faculty seminar. I covered my own work on the rudimentary framework for Benchmark Solutions pricing service (now BMRK on Bloomberg) but did not go into details on all the other efforts of everyone involved. Some of the material has already been covered in my blog post on Kalman sensitivities via the adjoint state method.

Monday, December 14, 2015

Domino dancing. A medium data research workflow with Ruffus on AWS.

I don't normally blog about technology because frankly, it is too much like work. But I have been trying Domino data labs of late, a nice product by a team mostly out of Bridgewater who clearly have good experience trying to make life easier for quants. They have built a research platform most readily used with Amazon compute but with the option of an in-house install. I certainly feel their efforts are worthy of note.

I don't really know what a data science platform should be necessarily, but perhaps it is a little like picking a dance partner: they have to know when to get out of the way in order for things to move effortlessly. That was a quality I found lacking in both in-house and competing platforms that seem to assume you are trying to step in a certain direction and in the process, inadvertently make it difficult to do simple things you feel you have already mastered. In short, they feel a bit like dancing with ... err ... me.

My workflow

I'll describe my setup for context, though I doubt there is anything atypical about my needs. I had already decided on Amazon over some in-house technology, incidentally, for several key reasons including the on-demand capability and availability of an attached file system. Another preference was use of a single beefy machine rather than many small ones because some of my steps require all the data together. With 8 GPU machines standard, and a 100 processor option coming in 2016 it is hard to see the merit in administering my own box.

My econometric/algo work (a.k.a. data science these days) pushes thousands of data files through a few dozen processing steps. After running local tests I sync both data and code using the Domino client to AWS, then kick off a job there. When the job is done, I simply sync again and examine the results locally.

I use the free python library Ruffus to create a non-blocking pipeline and handle the splitting and recombining of data files therein. I can't show you my actual computation graph as it is proprietary, so here's a stock image from ruffus instead:

Imagine the same with a few dozen steps and you get the idea. Much of what I do lives on one or more of these growing DAGs.

Serving results behind a REST API

The nodes on my computation tree comprise between one and a few thousand files of the same type. The leaves comprise summary tables (sometimes documentation) for this and that. Some output becomes lookup tables for REST API's deployed on Amazon. Domino provides simple version management for the code and data serving those endpoints, making this effortless.

One merely provides the entry point: in my case a plain python function to respond to the request (no wrapping/unwrapping). At present Domino only allows a single entry point but that is hardly a restriction. I merely register the following dispatcher function to automatically exposes every other function I drop in the same module.

That's REST, pretty much for free.

Monitoring runs

No complaints here. Domino presents your runs to colleagues alongside run diagnostics such as the following.

It could be slightly easier to arrange the results, arguably, but one way around that is to maintain your own markdown index page ( and point users to where you really want them to be. As a minor side point I discovered that the markdown rendering was just slightly non-standard - but I'm sure that will be knocked off the issue list by the time you read this.

Serving richer results

By default Domino tracks the new files produced by your runs and renders them to collaborators browsing the results. If you want to serve up richer content than a pile of .csv's you can have your run create an .html file and use javascript. As suggested by one of the Domino developers, here is one way of adding sortability to output tables:

Thanks John. The same goes for plots. Need something fancy? Just add the package to your requirements.txt and it will be installed prior to the run starting.

Reproducible research

In a sense Domino is pioneering what truly reproducible research might look like. I made a very poor stab at the same vague concept many years back so I am sympathetic as to the difficulties. (That effort lived on this site incidentally, and was so dismal I'd forgotten it myself until just now.)

There is an opportunity. While the workflow for developing code is obviously well established in source control programs and accompanying verbiage, the organizing of code, data and intermediate temporary data tends to be more bespoke. The Domino starting point, philosophically, is that code and data be kept in sync and this ensures absolute reproducibility of results. But they have not baked any religious assumptions into the product and one is free to break that symmetry. I'll get to why it helped me to do so a little later.

Certainly the major annoyances (code sharing with package management) and organizational features (basic comment tracking etc) are cleanly solved if your colleagues are also using Domino. One points them to a link to start the run that reproduces your work, or they can select from a menu on your recent runs.

Transparent research

While multi-directional collaboration is a marketing point for Domino, I find that the weaker variety of this (call it transparency) is a sufficient motivation. There is certainly no barrier here for non-technical end-users who want to re-run your research with different parameters and then compare the results. And that comparison is straightforward. One selects two finished runs in the run history and selects the "compare" checkbox to get the full diff.

Note however, that if you keep your code in the same project this comparison will include all your source changes. If you don't want that for yourself or your readers, a simple hack here is to use .dominoresults to hide your code and echo major meta-parameters to their own output file in lieu. Then the Domino presentation of your backtest results (say) is actually very helpful.

I'm reasonably happy with the level of transparency that the Domino platform provides. I feel it will replace a fair amount of presentation preparation (I'm a bit lazy on that front). I maintain a stable version of data and code for this purpose - or at least I aspire to - so that others can drill down into very large runs. One might hope that the days of manually converting piles of data into long powerpoint presentations are behind me, though that might be optimistic.

Perhaps an area for improvement or enhancement of the Domino platform would involve documentation. Some support for Latex generation of documents for instance.

Separating code and data again

As my data sets grew, I ultimately found it more convenient to use a separate project for code and several different projects for data. There are distinct categories of data living remotely and locally. Some data needs to be synced between local and remote instances to make it easier to run quick experiments. Some data is merely temporary. Some data should only ever live on Amazon due to it's size.

Then there is the pragmatic concern with code. I personally didn't want the possibility of large globs of data corrupting my defacto code repo. Nor was there any need for the bulk of my code to be shared with those wanting to make minor hacks to the results.

These concerns are addressed. The Domino team was very helpful and introducing read-only sharing between projects. Thereafter I moved to a pattern of "code only" and "code-light data heavy" projects. In my project holding source code I allow read-only export to the other projects.

To minimize syncing of large files (and lots of files) to amazon I use a chain of Domino data projects. My medium-data flow is as follows:

  1. Ingest data locally to ~/zipped_data_project/
  2. Sync the compressed data to Amazon with Domino client
  3. Unzip on Amazon to a /mnt//unzipped_data_project/... by allowing this second project to import files from /zipped_data_project, (and also from the source project containing the script in question and the rest of one's code).
  4. Allow a results project /mnt//endpoint_project to read the unzipped data and run the pipeline, creating result files and lookup tables. Use .dominoignore config to decide what is worth keeping.
  5. Publish the endpoint.
I've leaving out source control & backup since there are many ways to do that. I use my own hacks to ensure I always backup right before syncing, because very occasionally one may need to recover from a git issue under the hood. The implicit version control achieved by syncing with the Domino client is, for me, more of a secondary line of defense.

I suspect some future improvements to Domino might obviate the minor inconvenience of the zip/unzip shuffle. And some of my issues there relate mostly to sub-par connectivity.

Overall impression

Sweet. Hopefully not sweet as in the Pepsi versus Coke trick (i.e. ask me again when I've drunk the whole can) but even if it is regarded as mere sugar for Amazon compute, Domino data labs seems to me to represent good, very reasonably priced sugar that doesn't get in the way of itself or you. It isn't trying to replace anything, just bring together a few critical ingredients like git, docker, jupyter, python/R/Julia package management and some friends from the Apache family.

The Domino team has cut the time for basics like REST APIs, job management, launchers for non-technical users, scheduling, sensible run management and results rendering. It is drawn together seamlessly with Scala (okay the fact I know that is a tell. But let's say 99% seamless with only the occasional stack trace finding its way to the user :-). You don't need a separate AWS account and overall, it is very simple.

Many quants or supporting groups might no doubt build their own versions of the same but I've come to realize that pulling all this together, or even a poor man's version, can easily fall foul of Hofstadter's Law. It is a little too soon to say that Domino have commoditized this area, but they seem to be off to a good start and with a decent size engineering team in place, I expect Domino to get better. As a quant, that's good news. Liberating even.

Costs seem to be well below internal costs. And my recommendation is that even for hobby projects involving only python, R and/or Julia one might want to consider whether pricing like the following is unreasonable when it obviates a great deal of messing around.

That's my 2 cents for now (or my 9.3 cents as the case may be). For those who are interested I've written up a few more hacks in this followup post.

Wednesday, October 29, 2014

Some special cases of the Lambert W function

Suppose \(f(x)\) is given by the functional relation \begin{equation} e^{-a(x) f(x) + b(x)} = c(x) f(x) + d(x) \end{equation} then \begin{equation} f(x) = \frac{1}{a(x)}W\left( \frac{a(x)}{c(x)} e^{b(x)+\frac{a(x)d(x)}{c(x)} } \right) - \frac{d(x)}{c(x)} \end{equation} where \(W\) is the Lambert W Function. You can click through to Wikipedia for a proof. Here I merely list some special cases which every good citizen should instantly recognize (and it can be mildly tedious to reproduce them).

Recovering the definition of \(W(x)\)

Just to check we are on the right planet set \(a=1\), \(c=1\), \(b=\ln(x)\) and \(d=0\) to get \begin{equation} x e^{-f(x)} = f(x) \end{equation} or equivalently \( f(x)e^{f(x)} = x\) which is the definition of the Lambert W function. Sure enough those substitutions recover $$ f(x) = W(x) $$ as expected. A minor modification uses \(a(x)=1\), \(c(x)=1\), \(b(x)=0\) but leaves \(d(x)\) general.

The solution to \( e^{-f(x)} - f(x) = d(x) \)

is \begin{equation} f(x) = W \left( e^{d(x)} \right) \end{equation} Only slightly more generally, set \(a(x)=a\), \(b=\ln(g(x)), c=1, d=1\) for those cases where \(g(x)\) is on the wrong side, as it were.

The solution to \( e^{-af(x)} = f(x) g(x) \)

or equivalently \( f(x) e^{af(x)} g(x) = 1 \) is \begin{equation} f(x) = \frac{1}{a} W\left( a g(x) \right) \end{equation} In particular,

The solution to \( x^k e^{-af(x)} f(x) = 1 \)

which seems to crop up a fair bit for your author is \begin{equation} f(x) = \frac{1}{a} W\left( a x^{-k} \right) \end{equation}

Similarly setting \(a=-1\) ...

The solution to \( x e^{f(x)} f(x) = 1 \)

(i.e. where \(x\) is on the wrong side but we otherwise have the Lambert W definition) must be \begin{equation} f(x) = - W\left( -\frac{1}{x} \right) \end{equation} We might also take \(b = \ln(g(x))\) and thus

The solution to \( g(x) e^{-af(x)} = f(x) \)

or equivalently \( f(x)e^{af(x)} = g(x) \) is \begin{equation} f(x) = \frac{1}{a} W\left( a g(x) \right) \end{equation} which reduces to the Lambert W function for \(a=1\), as we expect. It is also pretty obvious from first principles, because if we multiply both sides by \(a\) we have \begin{equation} (af) e^{(af)} = ag \end{equation} and thus \(af = W(ag)\). Next suppose we want a power of \(f\) to appear. Let \(b=\beta(x)/k\), \(c = \gamma^{1/k}\), \(a = k/\alpha\) and \(d=0\). And then raise both sides to the power \(k\). It follows that...

The solution to \( e^{\alpha f(x) + \beta(x)} = \gamma f(x)^k \)

is \begin{equation} f(x) = \frac{\alpha}{k}W\left( \frac{k e^{\frac{1}{k}\beta(x)} }{\alpha \gamma^{1/k}} \right) \end{equation} and if, in particular, \(\beta(x) = -\ln(g(x))\) then

The solution to \( e^{\alpha f(x)} = g(x) \gamma f(x)^k \)

is \begin{equation} f(x) = \frac{\alpha}{k}W\left( \frac{k }{\alpha g(x)^{1/k} \gamma^{1/k}} \right) \end{equation} and if we take \(c = \frac{1}{\gamma}\) and \(g(x)=x\) and \(\alpha = \frac{s}{2}\) and \(k=2\) then p>

The solution to \( e^{-\frac{s}{2}f(x)} x f(x)^2 = c \)

is \begin{equation} f(x) = \frac{s}{2k}W\left( \frac{2k}{s} \sqrt{ \frac{c}{x}} \right) \end{equation}

Wednesday, August 13, 2014

A new use for the number 59575553 (ISDA shorten's 20 character Legal Entity Identifiers to 10 character Unique Trade Identifier Prefixes)

The International Swaps and Dealers Association (ISDA) faced a minor problem recently. In shortening their 20 character Legal Entity Identifiers (LEI's) into 10 character Universal Trade Identifier (UTI) prefixes they ran into collisions.

Somewhat randomly this came across my desk. So under proposal is an improvement for hashing LEI's into UTI prefixes: the modulus operation lifted from integers to the space of case insensitive alpha-numeric strings.

You can grab the python code here. Of course it's hardly a new idea.

Wednesday, May 28, 2014

Sensitivities of a Kalman Filter estimate with respect to all past observations

Consider a sequence of observations \(t_1,...,t_k\) at which a latent vector process \(x\) is observed indirectly, via an observation equation \begin{equation} y_{t_i} = H_i x_{t_i} + \epsilon_i \end{equation} We assume \(\epsilon_i\) is mean zero multivariate gaussian with covariance \(R_i\). For brevity we refer to \(y_{t_i}\) as \(y_i\), \(x_{t_i}\) as \(x_i\) and so forth. We assume the evolution of \(x\) in between the times specified can be written \begin{equation} x_{i+1} = A_i x_i + u_i \end{equation} where \(u_i\) are also gaussian. In this linear gaussian system the recursive estimation of \(x_t\) is achieved by the well known Kalman filter, and the contemporaneous impact of the next observation \(y_{k+1}\) is also (it is merely proportional to the Kalman gain).

But less well appreciated is a related computation, the derivatives of the Kalman filter estimate with respect to a past observation \(y_i\). This note establishes how said computation can be achieved by making two observations. The first is a re-representation of the Kalman estimate in the form of a weighted least squares problem (not dissimilar to the Duncan Horn representation). The second observation is that sensitivities of any weighted least squares problem can be computed using an adjoint trick.

Step 1: The Kalman filter solution as a (particular) least squares problem

We shall set up a least squares problem involving the current state \(x_k\) only, and all the previous observations. We argue that the solution to this problem is identical to the Kalman filter. Since the current estimate \(\hat{y}_k\) is a simple linear function of the current state \(x_k\), this allows us to compute the derivative of the current estimate with respect to all previous observations.

In the Kalman filter we assume a gaussian prior on the initial state \(x_0\). This can introduce annoying special cases in what follows, but we can clean up the notation instead by introducing: \begin{eqnarray} y_{-1} & = & H_{-1} x_{-1} + \epsilon_{-1} \\ x_0 & = & A_{-1} x_{-1} + u_{-1} \end{eqnarray} provided \(H_{-1}\) and \(A_{-1}\) are identity matrices, \(\epsilon{-1}\) is identically zero, \(y_{-1}\) is set equal to the mean of our prior and \(u_0\) adopts its covariance. With the boundary conditions cleaned up in this fashion we can invert the dynamical equations, assuming only that \(A\)'s have left inverses \(A^{-1}\), as follows: \begin{equation} x_j = A^{-1}_{j}\left( x_{j+1} - u_j \right) \end{equation} and then re-arrange the observation equations so that the only value of \(x_i\) that appears is \(x_k\). \begin{eqnarray} y_k & = & H_k x_{k} + \epsilon_k \\ y_{k-1} & = & H_{k-1} x_{k-1} + \epsilon_{k-1} \\ & = & H_{k-1} \left( A^{-1}_{k-1}\left( x_{k} - u_{k-1} \right) \right) + \epsilon_{k-1} \\ & = & H_{k-1} A^{-1}_{k-1} x_{k} - H_{k-1} A^{-1}_{k-1} u_{k-1} + \epsilon_{k-1} \\ y_{k-2} & = & H_{k-2} x_{k-2} + \epsilon_{k-2} \\ & = & H_{k-2} \left( A^{-1}_{k-2}\left( x_{k-1} - u_{k-2} \right) \right) + \epsilon_{k-2} \\ & = & H_{k-2} A^{-1}_{k-2} x_{k-1} - H_{k-2} A^{-1}_{k-2} u_{k-2} + \epsilon_{k-2} \\ & = & H_{k-2} A^{-1}_{k-2} \left( A^{-1}_{k-1}\left( x_{k} - u_{k-1} \right) \right) - H_{k-2} A^{-1}_{k-2} u_{k-2} + \epsilon_{k-2} \\ & = & H_{k-2} A^{-1}_{k-2} A^{-1}_{k-1} x_{k} - H_{k-2} A^{-1}_{k-2} A^{-1}_{k-1} u_{k-1} - H_{k-2} A^{-1}_{k-2} u_{k-2} + \epsilon_{k-2} \\ & \dots & \end{eqnarray} from which it is apparent that if we write \(Y = (y_k, y_{k-1}, y_{k-2},...,y_{-1} ) \) then \begin{equation} Y = G x_{k} + \eta \end{equation} where \(G\) is the concatenation of the coefficients of \(x_k\) given above and \(\eta\) is the gaussian random variable equal to the sum of \(u_k\)'s and \(\epsilon_k\)'s (again, with coefficients as above, leading to a non-trivial covariance structure).

Step 2. (A review of the adjoint trick)

Suppose \(x\) solves \(Qx = b(y)\). The adjoint trick can be used to compute the derivative of \(g(x)\) w.r.t. y. In particular, if \(y\) is the observation and \(x\) the solution of a generalized least squares problem with error covariance \(R\) we can cast it in this form by writing: \begin{eqnarray} g(x)& = & H x Q & = & H^T R^{-1} H \\ b(y) & = & H^T R^{-1} y \end{eqnarray} Consider now \begin{equation} f(x,y) = 0 \end{equation} where \begin{equation} f(x,y) = Q x - b(y) \end{equation} We use derivatives of \begin{equation} \tilde{g} = g - \lambda^T f(x,y) \end{equation} with respect to \(y\) as a means of computing derivatives of \(g\) with respect to \(y\). Note that \begin{equation} \frac{\partial \tilde{g}}{\partial y} = \frac{\partial g}{\partial x}\frac{\partial x}{\partial y} - \lambda^T \left( \frac{ \partial f }{\partial x }\frac{\partial x}{\partial y} + \frac{\partial f}{\partial y} \right) \end{equation} and this will simplify if we choose \(\lambda\) judiciously as a solution of \begin{equation} \frac{\partial g}{\partial x } = \lambda^T \frac{\partial f}{\partial x} \end{equation} which we call the adjoint equation. For then \begin{eqnarray} \frac{\partial \tilde{g}}{\partial y} & = & \frac{\partial g}{\partial x}\frac{\partial x}{\partial y} - \lambda^T \left( \frac{ \partial f }{\partial x }\frac{\partial x}{\partial y} + \frac{\partial f}{\partial y} \right) \\ & = & -\lambda^T \frac{\partial f}{\partial y} \\ & = & \lambda^T \frac{\partial b}{\partial y} \end{eqnarray} Now specializing to \begin{equation} g(x) = H x \end{equation} and \(b(y)\) as above we can solve for this convenient choice of \(\lambda\) by writing \begin{eqnarray} H & = & \frac{\partial g}{\partial x} \\ & = & \lambda^T \frac{\partial f}{\partial x} \\ & = & \lambda^T Q \\ & = & \lambda^T H^T R^{-1} H \end{eqnarray} where the second equality is the adjoint equation. It should be clear from this how to compute derivatives of \(\tilde{g}\) with respect to \(y\), and thereby compute derivatives of \(g\) with respect to \(y\).

Wednesday, November 20, 2013

Keeping punters log-happy: Some properties of a "pristine parimutuel" market clearing mechanism

Is there such a thing as a perfect probabilistic paradise for punters? A place where participants' prior probabilistic prophecies are pulled apart and placed into a philosophically pristine parimutuel posterior, without painful departure from previously perceived pseudo-optimal particulars? I believe the platonic possibility is promising (though to provoke, in practice we puritans may be prevented from partaking in this powerful portal whilst we are paralyzed by the paranoid parlance of the present day).

                                               A pristine market

Well now that I'm out of spit, here are some characteristics of a market for mutually exclusive, exhaustive outcomes that I would like to see:
  1. Every participant is forced to invest in a manner that optimizes their long run wealth. They will allocate their wealth to maximize \(E^P[\log {\rm wealth}] \) with respect to their best estimate of real world probability \(P\).
  2. Every participant's estimate of probability, on which 1. is evidently based, must include both private information (from which they might derive a probability measure \(R\) say) and the market probabilities themselves (call them \(Q\)). In contrast, to simply equate \(P\) with \(R\) in 1. is a dreadful, though rather common, mistake.
  3. Criteria 2. is achieved without the need for participants to monitor and respond to price feedback. In contrast the usual mechanism for including market information, as we find in racetrack parimutuels, is to provide provisional estimates of \(Q\) in realtime so that punters can react.
The last requirement motivated my ponderings on this topic. In our perceived utopia participants are given inventive to provide their own private probabilities \(R\) well before "race time" (as it were). They can then put their feet up and relax safe in the knowledge that 1. will be satisfied automatically.

Such is non-trivial because even investors whose bets are tiny with respect to the overall market volume (and thus have negligible price impact) must quite rationally react to partial price discovery even though in "theory", they should simply bet proportionately irrespective of the odds on offer (a lovely accident mentioned below). I use scare quotes because unless the odds are truly known, \(Q\) should really enter the optimization via the back door: an update to \(P\) acknowledging the market's ability to aggregate information.

In this setup the market clearing mechanism does the heavy lifting normally performed iteratively and imperfectly. So long as there is a rule governing the manner in which "\(R + Q = P\)" (figuratively speaking) we can anticipate every player's estimate of \(P\) given \(Q\), and thereby solve simultaneously for Q and log-optimal allocations that presumably influence Q. Here we consider the simplest rule of that sort. We shall assume P is a convex combination of \(R\) and \(Q\) with fixed coefficients.

Also for simplicity we consider the case where the market is not subsidized (though that might be an interesting direction for generalization). Then linear pricing forces us to adopt the most straightforward parimutuel clearing mechanism once we known the allocations: divide all the money wagered amongst those choosing the winning horse, in proportion to their bet.

                                       Necessary optimality conditions

Suppose market participant \(k\) allocates all his wealth \(W^k\) across \(I\) mutually exclusive outcomes. Suppose his estimate of probability for the \(i\)'th state is given by \begin{equation} p^k_i = \eta^k r^k_i + (1- \eta^k) q_i \end{equation} where \(r^k\) is his best estimate using only private information, and \(q_i\) is the market implied probability arrived at by means to be discussed.

As noted this equation is a statement of a seemingly rational philosophy, independent of how the market operates. Investor \(k\) might have noticed in the past that his private information adds some explanatory power to the market, but he probably shouldn't ignore the market prices altogether in arriving at the best estimate of real world probability.

We shall further suppose, in what follows, that all \(K\) participants are rational in another sense. They wish to optimize the log of their posterior wealth. Now it is well appreciated that if \(p^k_i\) are considered fixed this log-optimality leads simply to proportional betting, but that is not the case here. Only \(r^k_i\) and \(\eta^k\)'s are fixed, and we shall attempt to construct allocations \(u = \{u^k\}\) and clearing prices \(q_i\) that overtly depend on the investments made by participants.

To that end let \(u^k\) denote the fraction of wealth investor \(k\) invests in the \(i\)'th state. Suppose that the market clears in parimutuel fashion, meaning that all participants receive the same price. The market probability for the \(i\)'th state must be \begin{equation} q_i = \frac{ \sum_k u^k_i W^k } { \sum_{k=1}^K W^k } = \sum_k u^k_i W^k \end{equation} since we might as well suppose w.l.o.g. that \( \sum_k W^k = 1 \).

The question then arises, is there a choice of \(\{u^k_i \}_{i,k}\) such that investment by each participant is log-optimal? Intuitively one would expect a market equilibrium provided \(\eta^k\)'s are strictly between \(0\) and \(1\).

The utility function for the \(k\)'th investor is \begin{eqnarray} U^k( u^k ) & = & E^p\left[ \log\left( \frac{u^k_i}{q_i} \right) \right] \\ & = & \sum_{i=1}^I p^k_i \log\left( \frac{u^k_i}{ \sum_k u^k_i W^k } \right) \\ & = & \sum_{i=1}^I ( \eta^k r^k_i + (1- \eta^k) q_i ) \log\left( \frac{u^k_i}{ \sum_k u^k_i W^k } \right) \end{eqnarray} and by definition of \(u^k_i\) the constraints are \(\sum_i u^k_i = 1 \). Or for every \(k\) we might write \( g(u^k)=0\) where \(g(u^k) = \sum_i u^k_i - 1\). This sets up the first order Lagrangian equations for \( \Lambda(u,\lambda) = U(u) - \lambda \cdot g( u ) \) where we collect the components \(k=1,..K\). As usual for these problems we set \( \nabla \Lambda = 0 \) because for optimality the derivative of the objective function must be proportional to the derivative of the constraint function. This leads to equations of the form \begin{eqnarray} 0 & = & p^k_i( q_i(u)) \left( \frac{1}{u^k_i} - \frac{W^k}{q_i(u)} \right) - \eta^k W^k \log \left( \frac{u_i}{q_i(u)} \right) - \lambda^k, \ \ \ {\rm and} \\ 0 & = & \sum_i u^k_i - 1 \end{eqnarray} relating the allocations \(u^k_i\). Notice there are \(KI+K\) equations and \(KI+K\) free variables, including both the \(u^k_i\)'s and the \(\lambda^k\) multipliers (the \(W^k\) are fixed parameters and we have \(q_i(u) = \sum_k u^k_i W^k\) as noted above). The solution to this system of non-linear equations defines a pristine parimutuel clearing mechanism.

                                  Comparison to proportional betting

In contrast if we imagine that \(q_i\) are not a function of allocations \(u^k_i\) but fixed, and further suppose \(\eta^k = 1\) then we return to the overly stylized world where participants don't take market prices into account, preferring to believe their own homework is a sufficient statistic. The optimality conditions are instead: \begin{eqnarray} 0 & = & \frac{p^r_i}{u^k_i} - \lambda^k, \ \ \ {\rm and} \\ 0 & = & \sum_i u^k_i - 1 \end{eqnarray} It is apparent from the first equation that \(u^k_i \propto r^k_i\) and then, from the second, that \(u^k_i = r^k_i\). We recover the remarkable accident referred to as proportional betting, where the price \(q_i\) does not enter the picture. This works perfectly well for blackjack where \(r^k = p^k\) but not most real world games where markets inevitably supplement one's private information.


While I have caged this discussion in market language, it should be apparent that we have derived an algorithm for combining multiple probabilistic forecasts into a single probability distribution \(Q = \{q_i\}\), for which there is a large literature. See the references to Ranjan and Gneiting for example.

I shall concede that what is philosophically inconsistent is my use of simple linear pooling to arrive at the individual's subjective probability estimates \(p^k = \{p^k_i\}_{i=1}^I\) based on their prior \(r^k\) and the final market price \(q_i\), given that we then derive a more sophisticated combination scheme for meshing between individuals. Why not use the "better" scheme to combine \(Q\) and \(R\)?

I suppose a flippant response is "why not put a picture of a boy holding a cereal box on the cereal box?". Let me think about a more satisfying answer and get back to you.