## Thursday, 29 April 2010

### Waiting for software and school children

## Thursday, 22 April 2010

### Plotting Lambert's W Function

Here you go:

I've been developing code for computing complex values of Lambert's W function, and I thought I'd share the Riemann-surface plots above, which are, from top to bottom, Re(W), Im(W) and abs(W) for the branches -1 (Blues), 0 (Greens) and 1 (Oranges). The route I've taken to get the plots has been a good learning experience for me.

## First approach

My initial idea was to write a mex wrapper to my algorithmic Fortran code, along the lines of what happens for the NAG Toolbox for MATLAB. I also wanted to see if I could provide a vectorized interface, as native MATLAB functions have.Whew, it's not much fun debugging mex files, is it?

The main thing I learned from this experience was to use

`mxCreateNumericMatrix`to create my function's return values: in this case the W value, the residual error and the fail flag.

/* prhs[2] is the array of input zs */ dims = mxGetDimensions(prhs[2]); plhs[0] = mxCreateNumericMatrix(dims[0], dims[1], mxDOUBLE_CLASS, mxCOMPLEX); if (nlhs >= 2) { plhs[1] = mxCreateNumericMatrix(dims[0], dims[1], mxDOUBLE_CLASS, mxREAL); if (nlhs >= 3) { plhs[2] = mxCreateNumericMatrix(dims[0], dims[1], mxINT_CLASS, mxREAL); } }Other approaches I had tried, such as populating local

`mxArray *`s which were only assigned back to the

`plhs`arrays at the end (and if required) caused all sorts of upset to MATLAB when it exited my mex routine.

With the working code, I managed to get some good plots in a MATLABy way. I felt a little discontent though, and then I remembered the advice of my old (and sadly departed) A-level statistics teacher, about how it's better to be able to solve a problem two ways than one.

## Second approach

I thought back to Mike Croucher's PyNAG tutorials on calling the NAG C Library from Python, and decided it might be fun to call my algorithmic Fortran code from Python and use matplotlib to do the plotting. While I was at it I decided I'd see if I could get Fortran 2003's C Interoperability to help me make the call to Fortran from Python more straightforward; in the 'unassisted' case, calling Fortran from Python is very similar to calling Fortran from C, with the usual headaches that entails.The primary hoop I had to jump through to allow C interoperability was to map Fortran character strings to arrays:

... USE, INTRINSIC :: iso_c_binding, ONLY : c_char ... CHARACTER (kind=c_char), DIMENSION (6), & INTENT (IN) :: in_str CHARACTER (kind=c_char), DIMENSION (10), & INTENT (OUT) :: out_str CHARACTER (len=SIZE(in_str)) :: local_in_str CHARACTER (len=SIZE(out_str)) :: local_out_str INTEGER :: i ... DO i = 1, SIZE(in_str) local_in_str(i:i) = in_str(i) END DO ... ! Here's where all the real stuff happens. Operate on ! local_in_str and local_out_str. ... DO i = 1, SIZE(out_str) out_str(i) = local_out_str(i:i) END DO ENDThen with a similarly-vectorized and Pythonic wrapper around my Fortran call I generated the plots above using

from mpl_toolkits.mplot3d import \ Axes3D from matplotlib import \ cm import matplotlib.pyplot as plt import numpy as np fig = plt.figure() ax = Axes3D(fig) n = 100 lower = -1 upper = -lower threshold = 0.001 X1, Y1 = np.meshgrid(np.linspace(lower, upper, n), np.linspace(lower, -threshold, n)) X2, Y2 = np.meshgrid(np.linspace(lower, upper, n), np.linspace(threshold, upper, n)) cmaps = {-1: cm.Blues, 0: cm.Greens, 1: cm.Oranges} # Loop over some branches for b in [-1, 0, 1]: w = nag_LambertW(X1 + 1j*Y1, b=b, fail=1)[0] ax.plot_surface(X1, Y1, w.imag, rstride=1, cstride=1, cmap=cmaps[b]) w = nag_LambertW(X2 + 1j*Y2, b=b, fail=1)[0] ax.plot_surface(X2, Y2, w.imag, rstride=1, cstride=1, cmap=cmaps[b]) plt.show()I had to do two separate subplots for each branch to avoid

`matplotlib`interpolating across the discontinuity at the branch cuts. Being a bit of a

`matplotlib`newb, I'm now wondering

- is there a better way of telling
`matplotlib`to avoid discontinuities? - is it possible to apply a colour map continuously across two subplots?
- doesn't
`plt.title('Title')`work with surface plots?

## Thursday, 15 April 2010

### New product launch = busy, but happy marketing department!

A new product launch (or release) is always a busy time for any marketing department. Months of planning and preparing new literature, presentations, collating benchmark data, new advertising campaigns, exhibition stands etc, and then comes the fun part - announcing to the outside world.

As with all marketing activities, there's a plethora of best practice guides in achieving a successful product launch. Our goal in the marketing department at NAG is to ensure that our news reaches every single user, potential user and specific media.

*How do we achieve this goal?*

The relatively recent emergence of new communication methods, blogging, microblogging and social media have transformed how we communicate with our audience and has made this launch different from those previous to it. What's so great about these new channels to market is that we can communicate and engage directly and instantaneously with those most important to us as an organisation. This is priceless.

Using these new methods in conjunction with more traditional marketing such as press releases and email blasts will serve to create some of the "buzz" that's needed to achieve early mindshare. However, even with the slickest of launch operations, if the product itself falls short of expectation and fails to impress then all the marketing effort is wasted.

Thankfully at NAG, our products reflect our reputation for quality and correctness. All new numerical routines have to undergo stringent quality checks. This vital testing process results in a portfolio of world-class products that simply wouldn't reach launch if they weren't something we could be proud of.

The NAG Library for SMP & multicore launched on 15th April 2010.

## Friday, 9 April 2010

### Scratch

How early can you be taught to program? A thought that has struck me since my wife gave birth to our first born son a month ago. There is no denying that the thought processes that have to be learnt in order to program are invaluable and transferrable abilities. Skills such as logical thinking, reasoning, planning, modularity, that are applicable in all forms of jobs not just ones involving programming.

Having only learnt to program three years ago, I wish I had had more of a chance to program in my early years. So if ever my son is interested in programming I might start him off with a kids programming language called scratch that I heard about in a magazine.

Scratch is a pictorial drag and drop language that feels a bit like lego in style. Different types of variables are different shapes and only fit in the right shaped holes. It's not a very powerful language, but a good educational tool for the youngsters as it is easy to use, understand, and they can create a lot of fun interactive programs very easily.

I’m all for encouraging young ones to learn useful skills, and who knows, these new prodigy scratchers may well be updating and maintaining the NAG library in years to come.

## Thursday, 1 April 2010

### Easter egg or not?

While working on a demonstration of NAG random number generators in Excel and their application in option pricing using the Monte Carlo method I noticed a considerable discrepancy in results while using unscrambled Sobol quasirandom numbers. Well, to start with, how do you price a European option when the underlying asset follows geometric Brownian motion? Basically for each random number you simulate what the price of the asset will be in, say T=1 year, moving in one step from the starting time point t=0 to the end point t=T, calculate the payoff, and the average value of all payoffs will be our desired option price (Monte Carlo methods for option pricing). I’m not going to get into the details of Monte Carlo method, that will be covered in a paper I'm currently working on. The idea is that this simulation works very well and gives good convergence to the result obtained using the close-form formula, no surprises here. The stock prices have log-normal distribution, as expected. This is illustrated in the figure below. This histogram groups the stock prices from 20000 simulations.

One could however want to price Asian options, which are path-dependent. The payoff of such options depends not on the final price of the asset (as in European options), but on the average price of the asset throughout the life of the option. We therefore need to generate a path of stock prices, so we’ll move from t=0 to t=T not in one step as previously, but in many steps, say 252 (this is the average number of trading days per year). In practice we generate 252 random numbers, generate a path of the asset price and for the final step we calculate the payoff. This is done for every Monte Carlo simulation, so as you’ve probably noticed we need 252 times more random numbers than in case of the European option. The quantity is not an issue, the Sobol sequence is long enough. What happens is that in case of such pathwise Monte Carlo simulation the distribution is not log-normal, it doesn’t actually resemble any distribution, as can be seen below.

Why is that? Why do the quasirandom numbers fail? The simulation using the same data, but pseudorandom numbers instead of quasirandom, does return stock prices that are log-normally distributed (below).

What is going on here? Well, I know the answer, it took me some time to realise, and it would have taken more time if not the help of my colleagues Jacques and Robert. Do you, dear reader of this blog, know why it is like this? Is this just an easter egg, or an interesting property of quasirandom numbers? Happy Easter!