Monday, November 3, 2008

GIS with Haskell 1

It's time to bite the bullet and do some GIS Haskelling.

My first project is to develop a simple map server in Haskell. Here are the ingredients:
  • PostgresSQL + PostGIS
  • Some data to put into the database. For this I sourced some Australian suburb boundaries.
  • A library for manipulation GIS geometry, GEOS. In particular this provides functions to parse WKT strings from the database into geometry structures.
  • Haskell CGI package
  • Haskell bindings to the GEOS library.
  • Extension to HaXML adding SVG combinators.
The definition of the map is specified using a set of combinators resulting in a DSL looking a lot like the MapFile format of MapServer.

The following gives a map showing the suburbs centred on the city of Melbourne.
ex0 = map (connection "host=localhost user=postgres password=postgres dbname=australia"
`u` size 700 700
`u` extents 144.897079467773 (-37.8575096130371) 0.16821284479996734 0.1410504416999956
`u` layer ( table "suburbs"
`u` geometry "the_geom"
`u` klass ( style ( outlinecolour 255 0 0 1
`u` colour 100 255 100 1))

The resulting SVG file when viewed looks like:
The components that make up the map definition:
  • connection - supplies the database connection parameters,
  • size - is the size of the SVG output to be generated.
  • extents - is the extents of the map in world coordinates.
  • layer - this defines a layer. The table property defines the database table to use, the geometry property defines the column name to use. The klass parameter defines the style to use for drawing the geometry.
Layers are the key to this. Basic layers associate a geometry from the database with a style to be used when drawing the geometry.

A slightly more complex example with two layers is the following:
ex1 = map (connection "host=localhost user=postgres password=postgres dbname=australia"
`u` size 700 700
`u` extents 144.897079467773 (-37.8575096130371) 0.16821284479996734 0.1410504416999956
`u` layer ( table "suburbs"
`u` geometry "the_geom"
`u` labelitem "name_2006"
`u` klass ( style ( outlinecolour 255 0 0 1
`u` colour 100 255 100 1))
`u` label (colour 255 255 0 1))
`u` layer (table "suburbs"
`u` geometry "geomunion(the_geom)"
`u` klass ( style ( outlinecolour 0 0 0 1 `u` width 4)))

This is the same as before but with a new layer that provides a thick border around the edge of all the suburb boundaries:

Next steps are to source some population data, to colour code the suburbs depending on population and to include a legend.

Friday, September 26, 2008

Financial Contracts, Haskell and Probability

This article brings together the ideas presented in the paper 'How to write a financial contract' (HWFC) and Martin Erwig's FPF module.

We are going to deal with a simple but common situation in finance - if I have a contract where I am going to receive $100 dollars in 3 years time what is that 'contract' worth to me now. How much would I pay to obtain that contract? In order to calculate the worth we need to consider what else I would do with the money and the most obvious action is to deposit it into a bank account that attracts interest.

The question is reposed then as: if I put x into a bank account then what is x if the final amount after 3 years is $100. This is easy if the interest is fixed, not so easy if it varies.

This blogpiece will provide a fragment of the implementation of HWFC that answers the above.

As this is literate Haskell some preliminaries:

> module Main where
> import Probability

HWFC introduces the concept of a value process which is a function from time to a random variable. We shall equate a random variable with a probability distribution and a definition of a value process is:

> type PR a = Int -> Dist a

For our interest rate model let us say that from year to the next the interest rate can either stay the same, increase by 1% or decrease by 1% all with equal likelihood. We can express this as:

> interest :: Floating a => a -> PR a
> interest i n = (n *. one) i where one start = uniform [start+1/100,start,start-1/100]

The *. function allows us to repeat a random process n times. The process here is to start with an interest rate and to move to the next years rate.

If this year the rate is 10%, after a couple of years the distribution looks like:

interest 10 2
10.0 33.3%
9.99 22.2%
10.01 22.2%
9.98 11.1%
10.02 11.1%

Let us put that to one side and look at the contracts side of things. I will short circuit the approach in the paper and dive directly into the valuation

> data Obs a = O { evalObs :: PR a }
> konst k = O (\t -> certainly k)
> lift f (O pr) = O (\t -> fmap f (pr t))
> lift2 f (O pr1) (O pr2) = O (\t -> joinWith f (pr1 t) (pr2 t))
> date = O (\t -> certainly t)
> data Contract = C { evalContract :: PR Float }
> cconst k = C $ \ _ -> certainly k
> when o c = C $ disc (evalObs o) (evalContract c)
> at t = lift2 (==) date (konst t)
> zcb t x = when (at t) x
> whenFirstTrue :: PR Bool -> Int
> whenFirstTrue prb = f 0 where f i = if prb i == certainly True then i else f (i+1)
> baseRate = 10

This is a value process such that if when the first argument is true, return the second otherwise calculated the discounted value of the first argument.

> disc :: PR Bool -> PR Float -> PR Float
> disc prb prd t = if prb t == certainly True then prd t else let s = prd t
> t' = whenFirstTrue prb
> in discount baseRate s (t'-t)
> discount :: Floating a => a -> Dist a -> PR a
> discount int final time = let intspread = interest int time
> in joinWith (\i s -> s / (1+i/100)) intspread final

Lets start with a trivial example to make sure that things are working as planned

> ex1 = cconst 100

The value of this contract, as a random variable, is:

evalContract ex1 0

100.0 100.0

> ex2 = zcb 3 (cconst 100)

The value of this contact is:

evalContract ex2 0

90.90909 25.9%
90.900826 22.2%
90.91736 22.2%
90.89256 11.1%
90.92562 11.1%
90.88431 3.7%
90.93389 3.7%

The PFP library has a function to provide the expected value which can be ask of a distribution. The expected value of our contract is:

expected $ evalContract ex2 0


About Me

Melbourne, Australia
I work for GE in Melbourne Australia. All views do not necessarily represent GE.