Ivan Idris Blog

Numpy Strategies 0.1.2

So we are now in Sprint 12 of Project “NumPy Strategies”. Last week you finished the task of determining states and defining a Markov chain model. This took you ten minutes even though the estimate was 2 days. Hey, that’s the power of NumPy! The rest of the time you spent playing computer games/pool/pinball or reading your favorite book NumPy Beginner’s Guide.

Pop Quiz

Let’s see whether you remember what we decided the states should be:

• flat F, up U, bottom B, average A, rally R
• flat F, up U, down D
• None of the above.
• This is a trick question!

The next task is:

• Create the stochastic matrix with transition probabilities

We will estimate this task to 2 days again for obvious reasons.This task can be split into several subtasks, but these are of course not estimated.

1. Initialize the stochastic matrix to 0 values

We have 3 possible start states and 3 possible end states for each transition. For instance, if we start from a U state, we could switch to:

• U again
• F
• D

Initialize the stochastic matrix with the NumPy zeros function.

1

SM = numpy.zeros((3, 3))

2. Create a list containing the signs

This is just a normal Python list.

1

signs = [-1, 0, 1]

3. For each sign select the corresponding start state indices

Now the code becomes a bit messy. We will have to use actual loops! Maybe a better solution will occur to me one day. We will loop over the signs and select the start state indices corresponding to each sign. Select the indices with the NumPy where function.

1
2
3

for i in xrange(len(signs)):
   #we start the transition from the state with the specified sign
   start_indices = numpy.where(states[:-1] == signs[i])[0]

4. Continue to the next iteration if no start states are found

This is just straightforward Python code.

1
2
3
4
5

   N = len(start_indices)
 
   # skip since there are no transitions possible
   if N == 0:
      continue

5. Compute the end state values

The end state indices in our simple model can be found by just getting the next available position.

1
2

   #find the values of states at the end positions
   end_values = states[start_indices + 1]

6. Calculate the transition probabilities

We can now count the number of occurrences of each transition. Dividing by the total number of transitions for a given start state gives us the transition probabilities for our stochastic matrix. This is not the best method by the way, since it could be overfitting.

1
2
3
4
5
6

   for j in xrange(len(signs)):
      # number of occurrences of this transition 
      occurrences = len(end_values[end_values == signs[j]])
      SM[i][j] = occurrences/float(N)
 
print SM

The result for the stochastic matrix is:

1
2
3

[[ 0.47787611  0.          0.52212389]
 [ 0.          0.          0.        ]
 [ 0.42753623  0.          0.57246377]]

As you can see the transitions involving the flat no change state have 0 probability, while the other transitions have a probability close to 50%. This is I guess what we would have expected.

If you liked this post and are interested in NumPy check out NumPy Beginner’s Guide by yours truly. Thanks goes to David W. Lambert for his recent review on Amazon.com and contribution to the unofficial errata. As far as I know there are still 3 reviews in the pipeline for December. One of them was announced on The Glowing Python blog. This announcement also contains a review of the free sample chapter. I will be back next week with the steady state vector task.

Numpy Strategies 0.1.0

Happy New Year! So I was collecting evidence for my prediction of the gold bubble bursting early in 2011 and I found this article. Very short summary – 14 June 2011. What do I think? I think this is completeley incorrect. They are off by at least two days.

The plan for today is:

1. Define a volatility band for the log of stock prices and mean reverting strategy based on dips below the band.
2. Optimization for the holding period of shares.
3. Profit and champagne(optional).

Numpy 4 portfolio closed

The Numpy 4 portfolio was a nice little stable portfolio. I think that it had to do with the Jarque-Bera test prefiltering. Numpy 4 seemed to select mainly energy ETFs. I also did some experimentation with trailing stops. This also might account for its success.

Here is a summary of the closed Numpy portfolios until now.

Stops update

The Numpy 5 portfolio had some startup issues, but this week was pretty OK.

The band

The first step in defining the band is transforming to the logarithm of the stock price. Next, I find the local minima’s. For now, I define the local minimum as the point with the lowest value of the set containing this point and 5 earlier and 5 later points.

1
2
3
4
5
6
7
8
9
10

def locMins(arr, period=1):
   mins = []
 
   for i in range(period,len(arr) - period):
      themin = min(arr[i - period : i + period])
 
      if arr[i] <= themin:
         mins.append(i)
 
   return mins

Then, I fit a line through these minima. If the fit fails, Numpy returns an empty residuals array – so I check for that.

1
2
3
4
5
6
7
8
9
10
11
12

   logc = log(c[: int(argv[1]) ])
   mins = locMins(logc, 5)
 
   y = []
 
   for m in mins:
      y.append( logc[m] )
 
   (a,b,residuals) = fitline(mins, y)
 
   if len(residuals) == 0:
      continue

The band is formed by adding and subtracting the standard deviation of log of the close price with respect to the regression line. I also check whether a certain percentage of the historical prices falls within this band.

1
2
3
4
5
6
7
8
9

def within_band_ratio(arr, a, b, err):
   count = 0
 
   for i in range( len(arr) ):
      if arr[i] <= (a * i + b + err):
         if arr[i] >= (a * i + b - err):
            count += 1
 
   return count / float( len(arr) )

This is what the histogram of the “within band ratio” looks like. Seems to be skewed the right way, don’t you agree?

Now it’s time for the mean reverting strategy that selects dips below the band.

1
2
3
4
5
6
7
8
9
10

   STD = std(logc)
 
   last = log( c[ int(argv[1] ) + 1 ] )
   if last >= ( a * len(logc) + b - STD):
      continue
 
   wbr = within_band_ratio(logc, a, b, STD)
 
   if wbr < float( argv[2] ):
      continue

Optimal holding period

Common sense tells us that the optimal holding period would be less than a day. However, sometimes we don’t have enough data for that or whatever. The best choice is of course associated with the largest return, so I try to estimate this. O yeah, I did something with time averaging, ignore it for now.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

...
def expected_return( p, period ):
   indices = arange( 0, len(p), period )
   q = take( p, indices )
   r = get_returns( q )
 
   return geomean( r )
 
o,c = loadtxt(path( argv[ 1 ] + '.csv' ), delimiter=',', usecols=(3,6), unpack=True, converters={3: float_converter }) 
 
ers = []
twers = []
maxEr = expected_return( c, 1 ) 
ers.append( maxEr )
 
maxTwer = expected_return( c, 1 )
twers.append( maxTwer )
optPeriod = 1
optTWPeriod = 1
 
for i in range(2, 30):
   er = expected_return( c, i)
   twer = er / i
 
   ers.append( er )
   twers.append( twer )
 
   if er > maxEr:
      maxEr = er
      optPeriod = i
 
   if twer > maxTwer:
      maxTwer = twer
      optTWPeriod = i
...

Here are the plots of holding periods and expected returns.

This spreadsheet shows the optimal holding periods.

A Vogon Portfolio

Resistance is useless

Prostetnic Vogon Jetzt captain of the yellow constructor ship Grossebanana was feeling a bit sorry for himself. Jetzt denied his crew their shore leave, but that was not enough to cheer him up. Nobody appreciated his poetry. Not even the two hitchhikers, they had just caught.

See, see the sky
Marvel at its big pink depths.
Tell me, do you
Wonder why the goat ignores you?
Why its foobly stare
makes you feel sleepy.
I can tell you, it is
Worried by your shazbaz facial growth
That looks like
A milk.
What's more, it knows
Your fudge potting shed
Smells of frog.
Everything under the big sky
Asks why, why do you even bother?
You only charm cheeses.

It wasn’t that bad, was it now? The other problem was that the Prostetnic and the whole Vogon Construction Fleet would be out of work soon. You see, the Vogons were in the business of planet demolition. In the past tearing down planets was necessary to make space for hyperspace travel. But now, with the invention of a new type of spaceship, this was no longer needed. Jetzt was good at two things. Shouting and planet demolition. He had a brilliant career first as a space guard, Junior Shouting Officer, Just Shouting Officer, Senior Shouting Officer, Super Senior Shouting Officer and eventually captain of Grossebanana.

Jetzt did something unusual. He asked the hitchhikers for advice. First, they advised him jokingly to become a pirate, which provoked half hour of the Vogon shouting. When the captain calmed down, the freeloaders asked him whether he had considered trading. That was a fine profession in which shouting was encouraged.

Yes, the Prostetnic thought to himself, that was not such a bad idea. The Vogon’s victims encouraged by the sudden silence started talking about some kind of a trading strategy, that was discovered on an ancient planet, known throughout the Galaxy as being “mostly harmless”. Unfortunately, soon after the strategy became popular, the planet was demolished by the Vogons. A few hours of shouting later the hitchhikers and Jetzt created the following portfolio:

The captain pulled on the yellow fish stuck in his ear and shouted.

Resistance is useless!

THE END

Random links of general interest

• twill
• R-Forge TradeAnalytics
• R in Finance

If you liked this post and are interested in NumPy check out NumPy Beginner’s Guide by yours truly.