Category Archives: Julia

Pivot tables in DataFrames.jl 1.4

By: Blog by Bogumił Kamiński

Re-posted from: https://bkamins.github.io/julialang/2022/11/18/unstack.html

Introduction

Some time ago I have written a post about the unstack function
in DataFrames.jl. Since DataFrames.jl 1.4 release this function has received
a major update of offered functionality. Currently it is possible to
easily create pivot tables using unstack as I am going to show you today.

The example we will use is simulation analysis of graphs, as in May 2023
I am co-organizing
18th Workshop on Algorithms and Models for the Web Graph.
If you are interested in this topic please come and join us. If you would
like to present your work, here is the Call For Papers.

The post was written under Julia 1.8.2, Graphs.jl 1.7.4,
ProgressMeter.jl 1.7.2, Plots.jl 1.36.2, and DataFrames.jl 1.4.3.

Generating the data

We will consider Erdős–Rényi random graphs. These graphs take two
parameters n (the number of nodes in the graph) and p edge probability
(independently for each edge). It can be calculated thus, that the expected
number of edges in this graph is pn(n-1)/2, so expected average node degree
is d=p(n-1). So, given n and d we can calculate p=d/(n-1). In Graphs.jl
you can generate these graphs using the erdos_renyi function.

Connected components of a graph are equivalence classes of a
reachability relation between pairs of nodes. Informally: in a connected
component all nodes can be reached from each other and no nodes outside
from a connected component can be reached from it.

Let me explain it by example:

julia> using Graphs

julia> using Random

julia> Random.seed!(1234);

julia> er = erdos_renyi(7, 0.2)
{7, 4} undirected simple Int64 graph

julia> collect(edges(er))
4-element Vector{Graphs.SimpleGraphs.SimpleEdge{Int64}}:
 Edge 1 => 6
 Edge 3 => 7
 Edge 4 => 6
 Edge 5 => 6

julia> connected_components(er)
3-element Vector{Vector{Int64}}:
 [1, 4, 5, 6]
 [2]
 [3, 7]

We created a random graph on 7 nodes with edge probability equal to 0.2.
It has 4 edges. Note, that our graph is undirected, so although Graphs.jl
displays edges like this: 1 => 6, the edges are bi-directional.
We note that node 6 is connected to nodes 1, 4, and 5. So they
form one connected component. There is also an edge between nodes 3 and 7,
so this is a second connected component. Finally node 2 is isolated, giving
us a third component. We can easily find the connected components of a graph
using the connected_components function.

We will want to investigate how the size of the largest connected component
of the Erdős–Rényi random graph depends on n and d using simulation.

We generate the data for n having values [100, 1000, 10_000, 100_000]
and d having values 0.5:0.1:2.0. For each combination
of values we repeat the experiment 64 times.

First create an initial data frame with the setup of the experiment:

julia> using DataFrames

julia> df = allcombinations(DataFrame,
                            d=0.5:0.1:2.0,
                            n=[100, 1000, 10_000, 100_000],
                            rep=1:64)
4096×3 DataFrame
  Row │ d        n       rep   
      │ Float64  Int64   Int64 
──────┼────────────────────────
    1 │     0.5     100      1
    2 │     0.6     100      1
    3 │     0.7     100      1
    4 │     0.8     100      1
    5 │     0.9     100      1
  ⋮   │    ⋮       ⋮       ⋮
 4092 │     1.6  100000     64
 4093 │     1.7  100000     64
 4094 │     1.8  100000     64
 4095 │     1.9  100000     64
 4096 │     2.0  100000     64
              4086 rows omitted

Now we compute p for each experiment:

julia> df.p = df.d ./ (df.n .- 1);

julia> df
4096×4 DataFrame
  Row │ d        n       rep    p
      │ Float64  Int64   Int64  Float64     
──────┼─────────────────────────────────────
    1 │     0.5     100      1  0.00505051
    2 │     0.6     100      1  0.00606061
    3 │     0.7     100      1  0.00707071
    4 │     0.8     100      1  0.00808081
    5 │     0.9     100      1  0.00909091
  ⋮   │    ⋮       ⋮       ⋮         ⋮
 4092 │     1.6  100000     64  1.60002e-5
 4093 │     1.7  100000     64  1.70002e-5
 4094 │     1.8  100000     64  1.80002e-5
 4095 │     1.9  100000     64  1.90002e-5
 4096 │     2.0  100000     64  2.00002e-5
                           4086 rows omitted

Finally for each row we compute the fraction of nodes in the largest component
(the computation takes some time, so I use @showprogress to monitor it):

julia> using ProgressMeter

julia> df.largest = @showprogress map(eachrow(df)) do row
           er = erdos_renyi(row.n, row.p)
           cc = connected_components(er)
           return maximum(length, cc) / row.n
       end;
Progress: 100%|███████████████████████████████████████████████| Time: 0:00:39

julia> df
4096×5 DataFrame
  Row │ d        n       rep    p            largest 
      │ Float64  Int64   Int64  Float64      Float64 
──────┼──────────────────────────────────────────────
    1 │     0.5     100      1  0.00505051   0.04
    2 │     0.6     100      1  0.00606061   0.07
    3 │     0.7     100      1  0.00707071   0.22
    4 │     0.8     100      1  0.00808081   0.29
    5 │     0.9     100      1  0.00909091   0.35
  ⋮   │    ⋮       ⋮       ⋮         ⋮          ⋮
 4092 │     1.6  100000     64  1.60002e-5   0.64289
 4093 │     1.7  100000     64  1.70002e-5   0.69333
 4094 │     1.8  100000     64  1.80002e-5   0.73421
 4095 │     1.9  100000     64  1.90002e-5   0.77021
 4096 │     2.0  100000     64  2.00002e-5   0.79805
                                    4086 rows omitted

Now we are ready to analyze this data using unstack.

Analyzing the results

First let me recall the general structure of the unstack call
(please refer to its documentation for more details):

unstack([data_frame],
        [columns to use as row keys],
        [column storing column names],
        [column storing the values];
        combine=[function to apply to the values])

The structure is easy to remember. After a data frame, we sequentially pass
what we want in rows, what we want in columns, what values we want to unstack,
and finally, a combine keyword argument (introduced in DataFrames.jl 1.4)
specifies what operation we want to apply to these values.

By default combine is only, which means that we want a unique value in
row-column combination and otherwise we get an error. Let us check how
unstack works by default by using [:d, :rep] for rows, :n for columns,
and :largest for values:

julia> unstack(df, [:d, :rep], :n, :largest)
1024×6 DataFrame
  Row │ d        rep    100       1000      10000     100000   
      │ Float64  Int64  Float64?  Float64?  Float64?  Float64? 
──────┼────────────────────────────────────────────────────────
    1 │     0.5      1      0.04     0.009    0.0017   0.00021
    2 │     0.6      1      0.07     0.046    0.0024   0.00035
    3 │     0.7      1      0.22     0.011    0.0042   0.00048
    4 │     0.8      1      0.29     0.019    0.01     0.00139
    5 │     0.9      1      0.35     0.027    0.0109   0.00166
  ⋮   │    ⋮       ⋮       ⋮         ⋮         ⋮         ⋮
 1020 │     1.6     64      0.63     0.665    0.6419   0.64289
 1021 │     1.7     64      0.69     0.683    0.6941   0.69333
 1022 │     1.8     64      0.33     0.746    0.7389   0.73421
 1023 │     1.9     64      0.84     0.761    0.7696   0.77021
 1024 │     2.0     64      0.83     0.818    0.7991   0.79805
                                               1104 rows omitted

Everything worked, because for each cell we had a unique combination of row and
column keys. However, if we, for example, dropped :rep, we would get an
error:

julia> unstack(df, :d, :n, :largest)
ERROR: ArgumentError: Duplicate entries in unstack
at row 65 for key (0.5,) and variable 100.
Pass `combine` keyword argumen t to specify how they should be handled.

Let us start with passing combine=copy to just store all values of
:largest per row-column combination:

julia> unstack(df, :d, :n, :largest, combine=copy)
16×5 DataFrame
 Row │ d        100                                1000                              
     │ Float64  Array…?                            Array…?                          
─────┼───────────────────────────────────────────────────────────────────
   1 │     0.5  [0.04, 0.06, 0.05, 0.05, 0.08, 0…  [0.009, 0.015, 0.009,
   2 │     0.6  [0.07, 0.07, 0.05, 0.07, 0.04, 0…  [0.046, 0.013, 0.018,
   3 │     0.7  [0.22, 0.13, 0.05, 0.09, 0.18, 0…  [0.011, 0.019, 0.024,
   4 │     0.8  [0.29, 0.11, 0.08, 0.09, 0.1, 0.…  [0.019, 0.068, 0.018,
   5 │     0.9  [0.35, 0.19, 0.19, 0.19, 0.11, 0…  [0.027, 0.029, 0.057,
   6 │     1.0  [0.42, 0.12, 0.15, 0.22, 0.16, 0…  [0.11, 0.068, 0.042, 
   7 │     1.1  [0.26, 0.32, 0.23, 0.69, 0.36, 0…  [0.146, 0.085, 0.266,
   8 │     1.2  [0.45, 0.09, 0.35, 0.42, 0.25, 0…  [0.232, 0.175, 0.338,
   9 │     1.3  [0.49, 0.47, 0.12, 0.17, 0.67, 0…  [0.376, 0.504, 0.434,
  10 │     1.4  [0.36, 0.45, 0.49, 0.49, 0.47, 0…  [0.448, 0.527, 0.443,
  11 │     1.5  [0.72, 0.6, 0.65, 0.57, 0.69, 0.…  [0.577, 0.474, 0.581,
  12 │     1.6  [0.31, 0.77, 0.41, 0.42, 0.73, 0…  [0.639, 0.584, 0.644,
  13 │     1.7  [0.78, 0.62, 0.68, 0.7, 0.67, 0.…  [0.673, 0.651, 0.653,
  14 │     1.8  [0.86, 0.91, 0.7, 0.71, 0.76, 0.…  [0.721, 0.739, 0.748,
  15 │     1.9  [0.8, 0.76, 0.79, 0.83, 0.76, 0.…  [0.791, 0.794, 0.82, 
  16 │     2.0  [0.7, 0.74, 0.87, 0.88, 0.83, 0.…  [0.809, 0.776, 0.792,
                                                        3 columns omitted

Sometimes such operation is useful, but often we are interested in some
aggregates.

First let us check if indeed the data was generated correctly, i.e. if the
grid indeed has 64 experiments per combination of parameters. For this we use
the combine=length keyword argument (as we just want to check the number
of elements per dn combination:

julia> unstack(df, :d, :n, :largest, combine=length)
16×5 DataFrame
 Row │ d        100     1000    10000   100000 
     │ Float64  Int64?  Int64?  Int64?  Int64?
─────┼─────────────────────────────────────────
   1 │     0.5      64      64      64      64
   2 │     0.6      64      64      64      64
   3 │     0.7      64      64      64      64
   4 │     0.8      64      64      64      64
   5 │     0.9      64      64      64      64
   6 │     1.0      64      64      64      64
   7 │     1.1      64      64      64      64
   8 │     1.2      64      64      64      64
   9 │     1.3      64      64      64      64
  10 │     1.4      64      64      64      64
  11 │     1.5      64      64      64      64
  12 │     1.6      64      64      64      64
  13 │     1.7      64      64      64      64
  14 │     1.8      64      64      64      64
  15 │     1.9      64      64      64      64
  16 │     2.0      64      64      64      64

All looks good. So now we can check the mean and standard deviation of size of
the largest connected component for each combination:

julia> using Statistics

julia> unstack(df, :d, :n, :largest, combine=mean)
16×5 DataFrame
 Row │ d        100        1000       10000       100000      
     │ Float64  Float64?   Float64?   Float64?    Float64?
─────┼────────────────────────────────────────────────────────
   1 │     0.5  0.0525     0.0109375  0.00178125  0.000269219
   2 │     0.6  0.0776562  0.0159844  0.00266563  0.00040125
   3 │     0.7  0.0982813  0.0222969  0.00409688  0.000635625
   4 │     0.8  0.134375   0.031875   0.00668906  0.00121281
   5 │     0.9  0.161562   0.0485625  0.0136969   0.00312172
   6 │     1.0  0.205937   0.0905312  0.0445984   0.0187338
   7 │     1.1  0.265781   0.150547   0.149169    0.175325
   8 │     1.2  0.355625   0.272266   0.307183    0.313707
   9 │     1.3  0.375156   0.403047   0.421645    0.423478
  10 │     1.4  0.455313   0.503437   0.507891    0.510792
  11 │     1.5  0.594062   0.57925    0.580353    0.582628
  12 │     1.6  0.596719   0.63575    0.6428      0.641092
  13 │     1.7  0.697656   0.689906   0.691744    0.691409
  14 │     1.8  0.738281   0.730672   0.732986    0.732874
  15 │     1.9  0.76125    0.763125   0.768134    0.767024
  16 │     2.0  0.8025     0.797609   0.796216    0.796993

julia> agg_std = unstack(df, :d, :n, :largest, combine=std)
16×5 DataFrame
 Row │ d        100        1000        10000        100000      
     │ Float64  Float64?   Float64?    Float64?     Float64?
─────┼──────────────────────────────────────────────────────────
   1 │     0.5  0.0184305  0.00263598  0.000366829  4.57843e-5
   2 │     0.6  0.0422292  0.00599071  0.000623411  8.65108e-5
   3 │     0.7  0.0538238  0.00819539  0.00102446   0.000125444
   4 │     0.8  0.0635928  0.0138844   0.00234887   0.000429261
   5 │     0.9  0.0928426  0.0202656   0.00613093   0.00137704
   6 │     1.0  0.102581   0.0406045   0.026229     0.0111048
   7 │     1.1  0.141442   0.0748327   0.058437     0.0145289
   8 │     1.2  0.15247    0.0863053   0.0289068    0.00968446
   9 │     1.3  0.15323    0.0694691   0.0170002    0.00647732
  10 │     1.4  0.148399   0.0527765   0.0118846    0.00492798
  11 │     1.5  0.118987   0.0439715   0.0121344    0.00403464
  12 │     1.6  0.150407   0.0286661   0.00957977   0.00351468
  13 │     1.7  0.113833   0.0408923   0.00952214   0.00307108
  14 │     1.8  0.0969627  0.021692    0.00890106   0.00240749
  15 │     1.9  0.0829324  0.0224821   0.00788557   0.00210828
  16 │     2.0  0.0755929  0.0213895   0.0066356    0.00233577

What we can notice (I will not go into mathematics of these results; however,
if you find such computations interesting please be sure to join us during
WAW2023 conference):

  • as we increase average node degree the average size of the largest component
    increases;
  • there seems to be a sharp change for average degree above 1; it is especially
    visible as we increase number of nodes in the graph;
  • the results have largest standard deviation around average degree equal to 1;
    it drops for low and high values of average degree; also the standard
    deviation in general decreases as we increase number of nodes in the graph.

Let me additionally visualize the last observation:

julia> using Plots

julia> plot(agg_std.d, Matrix(agg_std[:, 2:end]),
            xlabel="average degree",
            ylabel="std of largest component relative size",
            labels=permutedims(names(agg_std, Not(:d))))

Here is a generated plot:

Standard deviation of largest component relative size

Conclusions

I hope you will find the new combine capabilities of unstack useful.

In DataFrames.jl 1.5 release we plan to add even more capabilities to unstack
(like allowing to use multiple columns to generate columns, or allowing to
process multiple value columns). I will post about it when it is done.

Trend Following with ETFs

By: Dean Markwick's Blog -- Julia

Re-posted from: https://dm13450.github.io/2022/11/18/Trend-Following-with-ETFs.html

Trend following is a rebranded name for momentum trading strategies. It looks at assets where the price has gone up and buying them because it believes the price will continue to rise and likewise for falling prices where it sells. I’ll use this post to show you how to build a basic trend-following strategy in Julia with free market data.


Enjoy these types of posts? Then you should sign up for my newsletter. It’s a short monthly recap of anything and everything I’ve found interesting recently plus
any posts I’ve written. So sign up and stay informed!






Trend following is one of the core quant strategies out there and countless pieces are written, podcasts made, and threads discussed all over Twitter about how to build and use a trend following strategy effectively. This blog post is my exploration of trend-following and uses ETFs to explore the ideas of momentum.

Over the past few months, we have been experiencing the first actual period of sustained volatility and general negative performance across all asset classes. Stocks no longer just go up. This is driven by inflation, the strengthening of the dollar, and the rise in global interest rates it feels like always buying equities 100% isn’t a foolproof plan and diversifying can help. This is where trend following comes in. It wants to try and provide positive returns whether the market is going up or down and, give diversification in regimes like we appear to be in today.

As ever I’ll be using AlpacaMarkets.jl for data and walking you through all the steps. My inspiration comes from the Cantab Capital (now part of GAM) where they had a beautiful blog post doing similar: Trend is not your only friend. Since moving over to the GAM website it no longer looks as good!

using AlpacaMarkets
using DataFrames, DataFramesMeta
using Dates
using Plots, PlotThemes
using RollingFunctions, Statistics
theme(:bright)

ETFs vs Futures in Trend Following

In the professional asset management world trend following is implemented using futures. They are typically cheaper to trade and easier to use leverage. For the average retail investor though, trading futures is a bit more of a headache as you need to roll them as they expire. So ETFs can be the more stress-free option that represents the same underlying.

More importantly, Alpaca Markets has data on ETFs for free whereas I think I would have to pay for the futures market data.

So let’s start by getting the data and preparing it for the strategy.

We will be using three ETFs that represent the different assets classes:

  • SPY for stocks
  • BND for bonds
  • GLD for gold

We expect the three of these ETFs to move in a somewhat independent manner to each other given their different durations, sensitivity to interest rates, and risk profiles.

The stock_bars function from AlpacaMarkets.jl returns the daily OHCL data that we will be working with. You also want to use the adjustment="all" flag so that dividends and stock splits are accounted for.

spy = stock_bars("SPY", "1Day"; startTime = now() - Year(10), limit = 10000, adjustment = "all")[1]
bnd = stock_bars("BND", "1Day"; startTime = now() - Year(10), limit = 10000, adjustment = "all")[1];
gld = stock_bars("GLD", "1Day"; startTime = now() - Year(10), limit = 10000, adjustment = "all")[1];

We do some basic cleaning, formatting the time into a DateTime and just pulling the columns we want, Open, Close Next Open.

function parse_date(t)
   Date(string(split(t, "T")[1]))
end

function clean(df, x) 
    df = @transform(df, :Date = parse_date.(:t), :Ticker = x, :NextOpen = [:o[2:end]; NaN])
   @select(df, :Date, :Ticker, :c, :o, :NextOpen)
end

spy = clean(spy, "SPY")
bnd = clean(bnd, "BND")
gld = clean(gld, "GLD");

Joining it all into one long data frame gives us the model data going forward

allPrices = vcat(spy, bnd, gld)
allPrices = sort(allPrices, :Date)
last(allPrices, 6)

6 rows × 5 columns

Date Ticker c o NextOpen
Date String Float64 Float64 Float64
1 2022-10-31 SPY 386.21 386.44 390.14
2 2022-10-31 BND 70.35 70.36 70.58
3 2022-10-31 GLD 151.91 152.16 153.82
4 2022-11-01 SPY 384.52 390.14 NaN
5 2022-11-01 BND 70.26 70.58 NaN
6 2022-11-01 GLD 153.46 153.82 NaN
plot(plot(spy.Date, spy.c, label = :none, title = "SPY"),
plot(bnd.Date, bnd.c, label = :none, title = "BND", color = "red"),
plot(gld.Date, gld.c, label = :none, title= "GLD", color = "green"), layout = (3,1))

ETF Closing Prices

The prices of each ETF are on different scales. BND is between 70-80, GLD in the 100’s and SPY is in the 300-400 range. The scale on which they move is also different. SPY has doubled since 2016, GLD 80% increase, and BND hasn’t done much. Therefore we need to normalise both of these factors to something comparable across all three. To achieve this we first calculate the log-returns of the close-to-close move of each ETF. We then calculate the rolling standard deviation of the price series to represent the volatility which is used to normalise the log returns

\[\hat{r} = 0.1 \frac{\ln p_t – \ln p_{ti}}{\sigma}.\]

allPrices = @transform(groupby(allPrices, :Ticker), 
                      :Return = [NaN; diff(log.(:c))],
                      :ReturnTC = [NaN; diff(log.(:NextOpen))]);

We also calculate the returns of using the NextOpen time series as a way to assess the transaction costs of the trend-following strategy, but more on that later.

runvar calculates the 256-day moving variance, which we take the square root of to get the running volatility. There the normalisation step is a simple multiplication and division.

allPrices = @transform(groupby(allPrices, :Ticker), :RunVol = sqrt.(runvar(:Return,  256)));
allPrices = @transform(groupby(allPrices, :Ticker), :rhat = :Return .* 0.1 ./ :RunVol);

Dropping any NaNs removes the data points before we had enough observations for the 256-day volatility calculation and calculating the cumulative sum of the returns gives us an equity curve.

allPricesClean = @subset(allPrices, .!isnan.(:rhat ))
allPricesClean = @transform(groupby(allPricesClean, :Ticker), :rhatC = cumsum(:rhat), :rc = cumsum(:Return));

To check this transformation has worked we aggregate across each ticker.

@combine(groupby(allPricesClean, :Ticker), :AvgReturn = mean(:Return), :AvgNormReturn = mean(:rhat),
                                           :StdReturn = std(:Return), :StdNormReturn = std(:rhat))

3 rows × 5 columns

Ticker AvgReturn AvgNormReturn StdReturn StdNormReturn
String Float64 Float64 Float64 Float64
1 SPY 0.000372633 0.00343358 0.0122161 0.115383
2 BND -9.4164e-5 -0.00223913 0.00349196 0.113982
3 GLD 0.000214564 0.00288888 0.0086931 0.101731

Summarising the average and standard deviation of the return and normalised return shows the standard deviation is now close to 0.1 as intended.

This is where leverage comes in because bonds have lower volatility than stocks, we have to borrow money to increase the volatility of our bond investment. For example, let’s borrow £100 with £10 of collateral and invest the £100 in BND. If BND moves up by 1% it is now worth £101, we sell and pay back our loan and we are left with £11 which is a 10% return on our original investment. So even though the price only moved by 1% the use of leverage has amplified our return by 10x.

Plotting the cumulative log returns shows how they are on similar scales now.

plot(allPricesClean.Date, allPricesClean.rhatC, 
     group = allPricesClean.Ticker, legend=:topleft, title = "Normalised Cumulative Returns")

ETF Log Returns

The Trend Following Signal

With our data in a good shape, we can move on to the actual signal construction. Just like the Cantab article, we will be using the 100-day moving average. Using the RollingFunctions.jl package again we just have to set the window length to 100 and it will do the hard work for us. The actual signal is whether this running average is positive or negative. If it is greater than zero we want to go long as the signal is saying buy. Likewise, if the rolling average is less than zero then we want to go short, the signal is saying sell. So this simply means taking the sign of the rolling average each day.

If we didn’t want to go short, we just want to know when to buy or sell to cover we can simplify it further by just using the signal when it is positive. We will call this the long only signal.

Using data frame manipulation we can calculate the signal per day for each ETF.

allPricesClean = @transform(groupby(allPricesClean, :Ticker), 
                            :Signal = sign.(runmean(:rhat, 100)), 
                            :SignalLO = runmean(:rhat, 100) .> 0);

Evaluating the Trend Following Strategy

We’ve got our signal that says when to go long and short for each ETF. We need to combine the return of each ETF per day to get out the Trend Following return. As we are using log returns this is as simple as summing across the ETFs multiplied by the signal on each day. We have three ETFs so need to weight each of the returns by 1/3 otherwise when comparing to the single ETFs we would have 3x as much capital invested in the trend following strategy.

portRes = @combine(groupby(allPricesClean, :Date), 
           :TotalReturn = sum((1/3)*(:Signal .* :rhat)), 
           :TotalReturnLO = sum((1/3)*(:SignalLO .* :rhat)),
           :TotalReturnTC = sum((1/3) * (:Signal .* :ReturnTC)),
           :TotalReturnUL = sum((1/3) * (:Signal .* :Return)));

Again, plotting the cumulative returns shows that this trend-following strategy (dark blue) is great. I’ve massively outperformed just being long the SPY ETF. Even when we remove the shorting element (Trend Following – LO, red) this has done well.

portRes = @transform(portRes, :TotalReturnC = cumsum(:TotalReturn), 
                              :TotalReturnLOC = cumsum(:TotalReturnLO), 
                              :TotalReturnTCC = cumsum(:TotalReturnTC),
                              :TotalReturnULC = cumsum(:TotalReturnUL))

plot(portRes.Date, portRes.TotalReturnC, label = "Trend Following", legendposition = :topleft, linewidth=3)
plot!(portRes.Date, portRes.TotalReturnLOC, label = "Trend Following - LO", legendposition = :topleft, linewidth=3)

plot!(allPricesClean.Date, allPricesClean.rhatC, group = allPricesClean.Ticker)

Trend following results

It’s up and to the right which is a good sign. By following the trends in the market we can profit when it is both going up and down. This is without any major sophistication on predicting the direction either. Simply using the average produces enough of a signal to be profitable.

Let’s focus on just what has happened this year

portRes2022 = @transform(@subset(portRes, :Date .>= Date("2022-01-01")), 
            :TotalReturnC = cumsum(:TotalReturn), 
            :TotalReturnLOC = cumsum(:TotalReturnLO),
            :TotalReturnULC = cumsum(:TotalReturnUL))

allPricesClean2022 = @subset(allPricesClean, :Date .>= Date("2022-01-01"))
allPricesClean2022 = @transform(groupby(allPricesClean2022, :Ticker), :rhatC = cumsum(:Return))

plot(portRes2022.Date, portRes2022.TotalReturnULC, label = "Trend Following", legendposition = :topleft, linewidth = 3)
plot!(portRes2022.Date, portRes2022.TotalReturnLOC, label = "Trend Following - LO", legendposition = :topleft, linewidth =3)
plot!(allPricesClean2022.Date, allPricesClean2022.rhatC, group = allPricesClean2022.Ticker)

2022 Trend Following Results

The long-only (red line) staying flat is indicating that we’ve been out of the market since July. The trend-following strategy has been helped by the fact that it is all one-way traffic in the markets at the minute. It’s just been heading lower and lower across all the asset classes since the summer.

Implementation Shortfall and Thinking Practically

Everything above looks like a decent trend-following approach. But how do we implement this, or simulate an implementation? By this, I mean obtaining the trades that should be made.

The mechanics of this strategy are simple:

  1. Calculate the close-to-close return of an ETF
  2. If it is above the 100-day moving average buy, if it’s below sell or go short.

The price we can trade at though is not the close price, markets are closed, and you can’t place any more orders! Instead, you will be buying at the open on the next day. So whilst our signal has triggered, our actual price is different from the theoretical price.

Sidenote, it is possible to trade after hours, and in the small size, you might be able to get close to the closing price.

Given this is retail and the sizes are so small, we are going to assume I get the next day’s open price. We can compare our purchase price to the model price. The difference between the two is called ‘Implementation Shortfall’ and measures how close you traded relative to the actual price you wanted to trade.

If we are buying higher than the model price (or selling lower) we are missing out on some of the moves and it’s going to eat the performance of the strategy.

To calculate this Implementation Shortfall (IS) we pull out the days where the signal changed, as this indicates a trade needs to be made.

allPricesClean = @transform(groupby(allPricesClean, :Ticker), :SigChange = [NaN; diff(:Signal)])

trades = @subset(allPricesClean[!, [:Date, :Ticker, :o, :c, :Signal, :NextOpen, :SigChange]], 
             :SigChange .!= 0);

Then by calculating the difference between the next open and closing price we have our estimation of the IS value.

@combine(groupby(trades, :Ticker), 
        :N = length(:Signal), 
        :IS = mean(:Signal .* 1e4 .* (:NextOpen .- :c) ./ :c))

3 rows × 3 columns

Ticker N IS
String Int64 Float64
1 SPY 40 14.7934
2 BND 69 -6.45196
3 GLD 61 10.4055

For both SPY and GLD we have lost 15bps and 10bps to this difference between the close and the next open. For BND though we actually would earn more than the close price, again, this is down to the one-way direction of bonds this year. But overall, implementation shortfall is a big drag on returns. We can plot the equity curve including this transaction cost.

plot(portRes.Date, portRes.TotalReturnULC, label = "Trend Following", legendposition = :topleft)
plot!(allPricesClean.Date, allPricesClean.rc, group = allPricesClean.Ticker)
plot!(portRes.Date, portRes.TotalReturnTCC, label = "Trend Following - TC", legendposition = :topleft, linewidth = 3)

Equity Curve with Transaction Costs

What’s the old saying:

“In theory there is no difference between theory and practice – in practice there is” (Yogi Berra)

Implementation shortfall is just one transaction cost. There are also actual physical costs to consider too:

  • Brokerage costs
  • Borrow fees
  • Capacity

Let’s start with the simple transaction costs. Each time you buy or sell you are probably paying some sort of small fee. At Interactive Brokers, it is $0.005 with a minimum of $1. So if you trade $100 of one of the ETFs you will be paying 1\% in feeds. You could get rid of this fee completely with a free broker like Robinhood which is one way to keep the costs down.

This model requires you to go short the different ETFs, so you will also need to pay borrow fees to whoever you borrow the ETF from. This is free for up to $100k in Interactive Brokers, so not a concern for the low-capacity retail trader. You can’t short sell on Robinhood directly, instead, you’ll have to use options or inverse ETFs which is just complicating matters and overall not going to bring costs down.

All these fees mean there is both a minimum amount and a maximum amount that you can implement this model with. If you have too small of a budget your returns will just be eaten up by transaction costs. If you have a large amount of money, the implementation shortfall hurts. This is what we mean by capacity.

Next Steps in Trend Following

The above implementation is the most basic possible trend following. It only uses three assets when there is a whole world of other ETFs out there. It has a simple signal (running average) and uses that to either allocate 100% long or 100% short. There are several ways in which you can go further to make this better.

  • Expand the asset universe.

Including more ETFs and ones that are uncorrelated to SPY, BND and GLD gives you a broader opportunity to go long or short. I would look at including some sort of commodity, real estate, and international equity component as the next ETFs.

  • A better trend signal.

The simple moving average is great as there are no moving parts, but nothing is stopping it from being expanded to a more sophisticated model. Including another moving average with a faster period, say 5 days, and then checking as to whether this faster average is higher or lower than the slow average is also commonly used.

  • Better asset allocation using the signal

The current signal says {-1, 1} and no in-between. This can lead to some volatile behaviour where the signal might be hovering around 0 and leading to going long and short quickly around multiple days. Instead, you should map the signal strength onto a target position with some sort of function. This could mean waiting for the signal to get strong enough before trading and then linearly adding to the position until you max out the allocation.

Conclusion

Trend following is a profitable strategy. We’ve shown that using a simple rolling average of the returns can produce a profitable signal. When scrutinising it a bit further we find that it is difficult to achieve these types of returns in practise. The overnight price movements in the ETFs means that we trade at a worse price. This combined with the general transaction costs of making the trades makes it a hard strategy to implement and to rely on another quote, there is no such thing as a free lunch, you have to put in a remarkable amount of effort to scale this kind of strategy up. We then highlight several ways you could take this research further. So you now have everything at your fingertips to explore Trend Following yourself. Let me know anything I’ve missed in the comments below.

Geospatial Data Science with Julia

By: Júlio Hoffimann

Re-posted from: https://juliohm.github.io/science/geodatascience/

Geospatial Data Science with Julia
presents a fresh approach to data science with
geospatial data and the Julia programming language.
It contains best practices for writting clean, readable
and performant code in geoscientific applications
involving sophisticated representations of the (sub)surface
of the Earth such as unstructured meshes made of 2D and
3D geometries.

By reading this book, you will:

  1. Get a broader perspective on geospatial data
  2. Learn advanced geostatistical algorithms
  3. Reproduce practical open source examples

Most importantly, you will learn a set of geospatial features
that is much richer than the
simple features
implemented in traditional geographic information systems (GIS).

If you would like to support the book,
consider starring it on GitHub:
https://github.com/JuliaEarth/geospatial-data-science-with-julia