By: Dean Markwick's Blog -- Julia

Re-posted from: https://dm13450.github.io/2022/08/05/FPL-Optimising.html

One of my talks for JuliaCon 2022 explored the use of JuMP to optimise a Fantasy Premier League (FPL) team. You can watch my presentation here: Optimising Fantasy Football with JuMP and this blog post is an accompaniment and extension to that talk. I’ve used FPL Review free expected points model and their tools to generate the team images, go check them out.

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!

Now last season was my first time playing this game. I started with an analytical approach but didn’t write the team optimising routines until later in the season, by which time it was too late to make too much of a difference. I finished at 353k, so not too bad for a first attempt, but quite a way off from that 100k “good player” milestone. I won’t be starting a YouTube channel for FPL anytime soon.

Still, a new season awaits, and with more knowledge, a hand-crafted optimiser, and some expected points, let’s see if I can do any better.

A Quick Overview of FPL

FPL is a fantasy football game where you need to choose a team of 15 players that consists of:

• 2 goalkeepers
• 5 defenders
• 5 midfielders
• 3 forwards

Then from these 15 players, you chose a team of 11 each week that must conform to:

• 1 goalkeeper
• Between 3 and 5 defenders
• Between 2 and 5 midfielders
• Between 1 and 3 forwards

You have a budget of £100 million and you can have at most 3 players from a given team. So no more than 3 Liverpool players etc.

You then score points based on how many goals a player scores, how many assists, and other ways. Each week you can transfer one player out of your squad of 15 for a new player.

That’s the long and short of it, you want to score the most points each week and be forwarding looking to ensure you are set for getting the most points.

A Quick Overview of JuMP

JuMP is an optimisation library for Julia. You write out your problem in the JuMP language, supply an optimiser and let it work its magic. For a detailed explanation of how you can solve the FPL problem in JuMP I recommend you watch my JuliaCon talk here:

But in short, we want to maximise the number of points based on the above constraints while sticking to the overall budget. The code is easy to interpret and there is just the odd bit of code massage to make it do what we want.

All my optimising functions are in the below file which is will be hosted on Github shortly so you can keep up to date with my tweaks.

include("team_optim_functions.jl")

FPL Review Expected Points

To start, we need some indication of each player’s ability. This is an expected points model and will take into account the player’s position, form, and overall ability to score FPL points. Rather than build my expected models I’m going to be using FPL Reviews numbers. They are a very popular site for this type of data and the amount of time I would have to invest to come up with a better model would be not worth the effort. Plus, I feel that the amount of variance in FPL points means that it’s a tough job anyway, it’s better to crowdsource the effort and use other results.

That being said, once you’ve set your team, there might be some edge in interpreting the statistics. But that’s a problem for another day.

FPL Review is nice enough to make their free model as a downloadable CSV so you can head there, download the file and pull it into Julia.

df = CSV.read("fplreview_1658563959", DataFrame)

To verify the numbers they have produced we can look and the total number of points each team is expected to score over the 5 game weeks they provide.

sort(@combine(groupby(df, :Team), 
       :TotalPoints_1 = sum(cols(Symbol("1_Pts"))),
       :TotalPoints_2 = sum(cols(Symbol("2_Pts"))),
       :TotalPoints_3 = sum(cols(Symbol("3_Pts"))),
       :TotalPoints_4 = sum(cols(Symbol("4_Pts"))),
       :TotalPoints_5 = sum(cols(Symbol("5_Pts"))) 
        ), :TotalPoints_5, rev=true)

So looks pretty sensible, Man City and Liverpool up the top, the newly promoted teams at the bottom. So looks like the FPL Review knows what they are doing.

With that done, let’s move on to optimising. I have to take the dataframe and prepare the inputs for my optimising functions.

expPoints1 = df[!, "1_Pts"]
expPoints2 = df[!, "2_Pts"]
expPoints3 = df[!, "3_Pts"]
expPoints4 = df[!, "4_Pts"]
expPoints5 = df[!, "5_Pts"]

cost = df.BV*10
position = df.Pos
team = df.Team

#currentSquad = rawData.Squad

posInt = recode(position, "M" => 3, "G" => 1, "F" => 4, "D" => 2)
df[!, "PosInt"] = posInt
df[!, "TotalExpPoints"] = expPoints1 + expPoints2 + expPoints3 + expPoints4 + expPoints5
teamDict = Dict(zip(sort(unique(team)), 1:20))
teamInt = get.([teamDict], team, NaN);

I have to multiply the buy values (BV) by 10 to get the values in the same units as my optimising code.

The Set and Forget Team

In this scenario, we add up all the expected points for the five game weeks and run the optimiser to select the highest scoring team over the 5 weeks. No transfers and we set the bench-weighting to 0.5.

# Best set and forget
modelF, resF = squad_selector(expPoints1 + expPoints2 + expPoints3 + expPoints4 + expPoints5, 
    cost, posInt, teamInt, 0.5, false)

It’s a pretty strong-looking team. Big at the back with all the premium defenders which is a slight danger as one conceded goal by either Liverpool or Man City could spell disaster for your rank. Plus no Salah is a bold move.

To add some human input, we can look at the other £5 million defenders to assess who to swap Walker with.

first(sort(@subset(df[!, [:Name, :Team, :Pos, :BV, :TotalExpPoints]], :BV .<= 5.0, :Pos .== "D", :Team .!= "Arsenal"), 
     :TotalExpPoints, rev=true), 5)

So Doherty or Digne seems like a decent shout. This just goes to show though that you can’t blindly follow the optimiser and you can add some alpha by tweaking as you see fit.

Update After Two Game Weeks

What about if we now allow transfers? We will optimise for the first two game weeks and then see how many transfers are needed afterward to maximise the number of points.

model, res1 = squad_selector(expPoints1 + expPoints2, cost, posInt, teamInt, 0.5)
currentSquad = zeros(nrow(df))
currentSquad[res1[1]["Squad"] .> 0.5] .= 1

res = Array{Dict{String, Any}}(undef, length(0:5))

expPoints = zeros(length(0:5))

for (i, t) in enumerate(0:5)
    model, res3 = transfer_test(expPoints3 + expPoints4 + expPoints5, cost, posInt, teamInt, 0.5, currentSquad, t, true)
    res[i] = res3[1]
    expPoints[i] = res3[1]["ExpPoints"]
end

Checking the expected points of the teams and adjusting for any transfers after the first two free ones gives us:

expPoints .- [0,0,0,1,2,3]*4

6-element Vector{Float64}:
 162.385
 164.295
 167.987
 165.767
 164.084
 161.726

So making two transfers improve our score by 5 points, so seems worth it. If we go beyond two transfers, then we will pay a 4 point penalty, so it seems worth

So Botman and Watkins are switched out for Gabriel and Toney. Again, not a bad-looking team, and making these transfers improves the expected points by 5.

Shortcomings

The FPL community can be split into two camps, those that think data help and those that think watching the games and the players help. So what are the major issues with these teams?

Firstly, Spurs have a glaring omission from any of the results. Given their strong finish to the season and high expectations coming into the season this is potentially a problem.

Things can change very quickly. After the first week, we will have some information on how different players are looking and by that time these teams could be very wrong with little flexibility to change them to adjust to the new information. I am reminded of last year where Luke Shaw was a hot pick in lots of initial teams and look how that turned out.

How off-meta these teams are. It’s hard to judge what the current template team is going to be at these early stages in the pre-season, but if you aren’t accounting for who other people will be owning you can find yourself being left behind all for the sake of being contrarian. For example, this team has put lots of money into goalkeepers when you could potentially spend that elsewhere.

Some of the players in the teams listed might not get that many minutes. Especially for the cheaper players, I could be selecting fringe players rather than the reliable starters for the lower teams. Again, similar to the last point, there is are ‘enablers’ that the wider community believes to be the most reliable at the lower price points.

And finally variance. FPL is a game of variance. Haaland is projected to score 7 points in his first match, which is the equivalent to playing the full 90 minutes and a goal/assist. He could quite easily only score 1 point after not starting and coming on for the last 10 minutes and you are then panicking about the future game weeks. Relying on these optimised teams can sometimes mean you forget about the variance and how easy it is for a player to not get close to the number of points they are predicted.

Conclusion and What Next

Overall using the optimiser helps reduce the manual process of working out if there is a better player at each price point. Instead, you can use it to inspire some teams and build on them from there adjusting accordingly. There are still some tweaks that I can build into the optimiser, making sure it doesn’t overload the defence with players from the same team and see if I can simulate week-by-week what the optimal transfer if there is one, should be.

I also want to try and make this a bit more interactive so I’m less reliant on the notebooks and have something more production-ready that other people can play with.

Also given we get a free wildcard over Christmas I can do a mid-season review and essentially start again! So check back here in a few months time.

By: Dean Markwick's Blog -- Julia

Re-posted from: https://dm13450.github.io/2022/07/08/Machine-Learning-Property-Loans.html

Estateguru is a website that lets you lend money to property
developers. They are relatively short loans of about a year with an
interest rate of 10%. Having a high-interest rate means that the loan
is more likely to default, so you will end up receiving none of your
money. But using the data they provide from their website, can we
build a machine learning model that will help us choose loans that
won’t go bad?

I turned this blog post into a conference talk that you can watch –
Machine Learning Property Loans for Fun and Profit | Dean Markwick | JuliaCon 2023

Enjoy these types of posts? Then sign up for my newsletter.

There is a variety of information available with each loan offer. You know what country it is in, what type of property, the interest rate, and the amount of money they are asking for relative to the total property value. All of these variables will be used to set the interest rate and the higher the interest rate, the more likely the loan will go bad. Earning more interest on a loan is a consequence of the higher risk. But what if the people at Estateguru set the interest rate wrong? Can we get a better model of predicting when a loan goes bad and use that to only invest in the ‘higher’ quality loans?

I’ve done something like this before as an interview task, trying to predict when a balance sheet loan might default.

Fire up your Julia notebook and let’s get predicting!

Environment

I’m running Julia 1.7 and have the latest versions of all these packages:

using DataFrames, CSV, DataFramesMeta
using Statistics, StatsBase, CategoricalArrays
using Plots, StatsPlots
using GLM, MLDataUtils

Getting the Data

First, we need to get the actual data. On Estateguru’s statistics page (https://estateguru.co/portal/statistics/?lang=en) they have a button to download the loan book in excel format. So do that and convert it into a CSV. This makes its machine readable for Julia. When using the CSV package to pull in the data we want to normalize the names (normalizenames=true) to remove the spaces from the column headers.

rawData = CSV.read("../../../Downloads/loan_book_06_05.csv",
     DataFrame;
     normalizenames=true);

We now want to find out the good, the bad, and the ugly loans. This is contained in the status column.

@combine(groupby(rawData, "Status"), :n=length(:Status))

We are interested in the bad loans:

• Late, in default, full recovered, partially recovered

and the good loans:

• Repaid

The open, fully invested, closed, and funded are currently ‘live’ so we
can’t do too much about them until we are happy with the model and can
forecast the probability of default. We will use the open loans to
predict a probability of default at the very end to give us some
investment ideas.

Even though you get money back in the fully/partially recovered case,
that’s a best-case scenario, so we should be pessimistic and assume
that you don’t get anything back.

When subsetting the bad loans, I learn from my previous mistakes in making something accidentally quadratic.

goodLoans = @subset(rawData, :Status .== "Repaid")
goodLoans[:, "BadLoan"] .= 0

badLoans = rawData[findall(in(["Late", "In Default", "Full Recovered", "Partially Recovered"]),
              rawData.Status), :]
badLoans[:, "BadLoan"] .= 1

openLoans = @subset(rawData, :Status .== "Open");

modelData = vcat(goodLoans, badLoans)
mean(modelData.BadLoan)

0.08508809626128062

So we have a ‘bad loan’ rate of 8.5%, and the probability of a random loan going bad is less than the average offered interest rate, so this looks like a win.

a = @combine(groupby(modelData, ["Currency", "BadLoan"]), 
         :N = length(:Currency),
         :TotalNotional = sum(:Funded_Total_Amount),
         :AverageInterestRate = mean(:Interest_Rate))
a.NotionalPct = a.TotalNotional / sum(a.TotalNotional)
a

But in terms of notional, it’s more like a 17% default rate. So quite high in the grand scheme of things. If I had invested in every property proportional to the amount offered, I would have lost 17 cents per EUR invested! So whilst the interest rate of 10% looks high, this just shows how your losses can depend on your investing strategy.

But, if we can predict what loans might default, we can avoid them, and collect the juicy return of the good loans. Let us throw some data science at the problem and see if we can make a decent prediction. But first up, let’s explore the data.

Exploring the Estateguru Data

All good modeling starts with looking at the distribution of different variables. Starting with the quantitive ones, we want to understand the variation of interest rates, the loan relative to the property value, the property value, and the total amount asked for.

These last two are heavy-tailed so we take the logarithm of the values.

plot(
  histogram(modelData[!, "Interest_Rate"], label=:none, title = "Interest Rate"),
  histogram(modelData[!, "LTV"], label=:none, title = "LTV"),
  histogram(log.(modelData[!, "Property_Value"]), label=:none, title = "Property Value"),
  histogram(log.(modelData[!, "Funded_Total_Amount"]), label=:none, title = "Funded Total Amount")
    )

The interest rate is around 10%-11% and the average LTV is around 60%
with sensible distributions. So need to transform anything there. The
property values and funded amounts had large tails and so had to be
log-transformed to be plugged into our model.

histogram((modelData[!, "Loan_Period"]), label = :none, title = "Loan Period")

The loan period is concentrated around 12 and 18 months, so we should divide the value by 12 to convert it into years.

Looking at the marginal distribution of the interest rate and LTV we can see that it is mainly concentrated around those 11% return and 60% LTV properties.

@df modelData marginalhist(:Interest_Rate, :LTV)

Moving onto the qualitative variables, we want to see what factors come up the most as a proportion of the total data.

function fracPlot(data::DataFrame, column::Symbol)
    sData = @combine(groupby(data, column), :frac=length(:Currency)/nrow(data), :n=length(:Currency))
    sort!(sData, :frac, rev = true)

    bar(sData[!, column], sData[!, "frac"],   title=String(column), label = :none, orientation=:h)
end

plot(
  fracPlot(modelData, :Country),
  fracPlot(modelData, :Schedule_Type),
  fracPlot(modelData, :Suretyship_existence),
  fracPlot(modelData, :Loan_Type)
)

60% of all the properties are in Estonia and most have a suretyship. For those that don’t know (like me before I googled it), a suretyship is an agreement with another party to pay the loan if the original debtor can’t. So like a guarantor when you are renting your uni flat.

fracPlot(modelData, :Property_Type)

Property type looks like it contains two pieces of information, we can split on the “-“ and add another column to our data.

modelData[!, "Property_Type_A"] .= "N/A"
modelData[!, "Property_Type_B"] .= "N/A"

for i in 1:nrow(modelData)

    pt = strip.(split(modelData[i, "Property_Type"], "-"))
    modelData[i, "Property_Type_A"] = pt[1]
    
    if length(pt) == 1
        modelData[i, "Property_Type_B"] = pt[1]
    else
        modelData[i, "Property_Type_B"] = join(pt[2:end], " ")
    end
end

sData = @combine(groupby(modelData, ["Property_Type_A", "Property_Type_B"]), 
                 :n = length(:Currency),
                 :DefaultRate = mean(:BadLoan))

That Summer Cottage is annoying, I’m going to change it into residential. We potentially could do the same for the “Multi Family/ Residential” types, but as there are 74 of them, with a sensible default rate, it can be revisited later.

modelData[findall(modelData[!, "Property_Type_A"] .== "Summer Cottage"),
                  "Property_Type_A"] .= "Residential"

All of our variables should help us predict whether a loan will be repaid back or not. There is enough overlap between the qualitative features and sensible distributions in the numerical variables for us to start building our model.

Machine Learning in MLJ.jl

In R I would use the caret package
to fit a variety of machine learning models. It provides a common interface
to many different packages so you can easily iterate through different model types to see which one works best with your data. In Julia, the MLJ.jl package does the same. By setting up the data structures correctly you can fit and evaluate different types of models on your data through one interface.

In my model process, I will be using a linear model, the xgboost package, a random forest model, and a k-nearest-neighbour model. This covers all the bases, linear, and tree-based models, so should give us a reasonable model at the end of the fitting process.

using MLJ

To begin with, I subset my original model data frame to just the columns needed and perform the same transformations.

modelData2 = modelData[:, ["BadLoan", "Interest_Rate", "Property_Value", "Funded_Total_Amount",
                           "LTV", "Loan_Period", 
                           "Country", "Schedule_Type", "Property_Type_A", "Property_Type_B",  
                           "Loan_Type", "Suretyship_existence"]]

modelData2[!, "log_Funded_Total_Amount"] = log.(modelData2[!, "Funded_Total_Amount"])
modelData2[!, "log_Property_Value"] = log.(modelData2[!, "Property_Value"])
modelData2[!, "LTV"] = modelData2[!, "LTV"] /100
modelData2[!, "Loan_Period"] = modelData2[!, "Loan_Period"] /12

modelData2 = select(modelData2, Not([:Funded_Total_Amount, :Property_Value]));

MLJ requires each column to have a correctly specified type. So I need to assign each factor column the Mulitclass type.

modelData2 = coerce(modelData2,
                      :BadLoan => Multiclass,
                      :Country=>Multiclass,
                      :Schedule_Type => Multiclass,
                      :Property_Type_A => Multiclass,
                      :Property_Type_B => Multiclass,
                      :Loan_Type => Multiclass,
                      :Suretyship_existence => Multiclass);

y, X = unpack(modelData2, ==(:BadLoan); rng=123);

train, test = partition(eachindex(y), 0.7, shuffle=true); # 70:30 split

The columns now have the correct types but need to now transform the multi-class X’s into integers, which is the equivalent of one-hot-encoding. This means training a ContinuousEncoder on the data and using it to expand the multi-class columns into multiple columns.

encoder = ContinuousEncoder()
encMach = machine(encoder, X) |> fit!
X_encoded = MLJ.transform(encMach, X)

We also want to scale the 3 numeric values to be mean 0 and standard deviation 1. This is the Standardizer model and again trained, but this time only on the training rows. We don’t want to leak information into the test set.

standardizer = @load Standardizer pkg=MLJModels
stanMach = fit!(machine(standardizer(features = [:Interest_Rate, :log_Funded_Total_Amount, :log_Property_Value]),
        X_encoded); rows=train)
X_trans = MLJ.transform(stanMach, X_encoded)

By using this machine-based workflow from MLJ we can make our transformations
repeatable and ensure no leakage from the train set to the test set.

With the data all prepared, we can now move on to fitting the models.

A Null Model

Like always, we need the null model to give our baseline performance. Our
predictions need to outperform this model to make sure the models are learning something.

With MLJ you need to pull in the model from an outside package, so you
will see this as a common pattern throughout this post.

constantModel = @load ConstantClassifier pkg=MLJModels

With it loaded we create a machine (that learns) with the data and then
evaluate its performance for several different measures.

We tell the machine to only use the train rows when fitting the
model. The resulting measures are the averages across the
cross-validation folds.

constMachine = machine(constantModel(), X_trans, y)

evaluate!(constMachine,
        rows=train,
         resampling=CV(shuffle=true),
         measures=[log_loss, accuracy, kappa, brier_loss, auc],
         verbosity=0)

We pass through some different metrics in the evaluation phase which
gives us some indication of how the model is performing. In this case,
a $\kappa$ of zero shows that the model isn’t doing anything more than
just predicting each loan will be fine, but we have a 91%
accuracy. This is because of the class imbalance, so we want to pay
attention to the Brier loss and area under the curve (AUC) to evaluate the model.

An Interest Rate Only Model

Next up is a model that just looks at the interest rate variable as a
predictor of the default. From some underlying maths, you can prove
that the default rate of a loan is proportional to the interest rate
offered. Higher interest rate loans have a higher risk, therefore they
have a higher reward, you need to be compensated for taking on this
higher probability of defaulting. We can fit a model that uses just
the interest rate column and see if we do any better than the null
model.

The beauty of MLJ means that we can just pull in the model type and
train a new machine just like previously.

This model is just a basic logistic classifier with two parameters,
the intercept, and the interest rate. We turn off the penalisation
(lambda=0) and run the model.

logisticClassifier = @load LogisticClassifier pkg=MLJLinearModels verbosity=0

irMachine = machine(logisticClassifier(lambda = 0),  X_trans[!, ["Interest_Rate"]], y)
fit!(irMachine, rows=train, verbosity=0)
evaluate!(irMachine, rows=train, 
          resampling=CV(shuffle=true), measures=[log_loss, accuracy, kappa, brier_loss, auc])

The AUC improves, but nothing else compared to the null model. So not
worth dwelling on really.

fitted_params(irMachine)

(classes = CategoricalValue{Int64, UInt32}[0, 1],
 coefs = [:Interest_Rate => 0.36875405640649356],
 intercept = -2.468541580016672,)

A positive value on the Interest_Rate coefficient confirms this. There
is an increase in the probability of default when the interest rate is higher.

A Linear Model

MLJ has one interface for both ridge/lasso regressive and logistic
regression. For logistic regression, we just set the penalisation value
(lambda) to 0 and train on the data but with all the features now.

lmMachine = machine(logisticClassifier(lambda=0), X_trans, y)

fit!(lmMachine, rows=train, verbosity=0)

evaluate!(lmMachine, 
          rows=train, 
          resampling=CV(shuffle=true), 
          measures=[log_loss, accuracy, kappa, brier_loss, auc],  verbosity = 0)

Now we are seeing some differences. The accuracy has increased and the kappa measure is now reporting a nonzero measure. So this model has learned something about the underlying data.

Penalised Regression

Now we can start to constrain the parameters of the linear model and
see if we can improve on the basic logistic regression. The
penalty=:en enables elastic-net regression, so both a \(L_1\) and
\(L_2\) penalisation that can help prevent the model overfitting.

lmModel = logisticClassifier(penalty=:en)

gamma_range = range(lmModel, :gamma, lower = 0, upper = 0.1)
lambda_range = range(lmModel, :lambda, lower = 0, upper = 0.1)

lmTuneModel = TunedModel(model=lmModel,
                          resampling = CV(nfolds=6, shuffle=true),
                          tuning = Grid(resolution=25),
                          range = [gamma_range, lambda_range],
                          measures=[auc, log_loss, accuracy, kappa, brier_loss]);

lmTunedMachine = machine(lmTuneModel, X_trans, y);

fit!(lmTunedMachine, rows=train, verbosity=0)

We can plot the results of the tuning of the different
hyperparameters.

plot(lmTunedMachine)

Not much performance differences over the different hyperparameters
which we can also see in the evaluation metrics, there is an actual
drop in performance. One such explanation could be the lack of overall
features compared to the number of observations, we don’t have a
massive amount of columns describing the loans.

So overall, this elastic-net model performs very little shrinkage and
doesn’t improve on the basic logistic model.

XGBoosting

We can now move on to tree-based models and the old faithful XGBoost. We will fit an untuned model and another model where we vary the hyperparameters.

xgboostModel = @load XGBoostClassifier pkg=XGBoost verbosity = 0

xgboostmodel = xgboostModel(eval_metric=:auc)

xgbMachine = machine(xgboostmodel, X_trans, y)

evaluate!(xgbMachine,
        rows=train,
         resampling=CV(nfolds = 6, shuffle=true),
         measures=[log_loss, accuracy, kappa, brier_loss, auc],
         verbosity=0)

Best model so far and without any tuning! But, we can see if we can vary
some of the hyperparameters and get a better fitting model.

I’ll vary \(\gamma, \eta, \lambda\) and \(\alpha\) from 0 to 5. Have
look at the
XGBoost docs
for an overview of what the parameters mean.

gamma_range = range(xgboostmodel, :gamma, lower = 0, upper = 5)
eta_range = range(xgboostmodel, :eta, lower = 0, upper = 1)
lambda_range = range(xgboostmodel, :lambda, lower = 1, upper = 5)
alpha_range = range(xgboostmodel, :alpha, lower = 0, upper = 5)

xgbTuneModel = TunedModel(model=xgboostmodel,
                          resampling = CV(nfolds=6, shuffle = true),
                          tuning = Grid(resolution=10),
                          range = [gamma_range, eta_range, lambda_range, alpha_range],
                          measures=[log_loss, accuracy, kappa, brier_loss, auc]);

xgbTunedMachine = machine(xgbTuneModel, X_trans, y);

fit!(xgbTunedMachine, rows=train, verbosity=0)

So a slight improvement in the \(\kappa\) metric, but the fact some of
the hyperparameters have gone to zero makes me think it’s going to be
overfitting the data. We will have to wait to look at the test
data performance.

K Nearest Neighbours

Next up I’ll use the K Nearest Neighbours algorithm to model the
data. This splits the data into \(K\) different chunks and attempts to
find the properties that are most similar and whether they default.

knnModel = @load KNNClassifier pkg=NearestNeighborModels verbosity = 0
knnmodel = knnModel()

knnMachine = machine(knnmodel, X_trans, y)
evaluate!(knnMachine,
        rows=train,
         resampling=CV(shuffle=true),
         measures=[log_loss, accuracy, kappa, brier_loss, auc],
         verbosity=0)

Five is the default number of neighbours, but that is a
hyper-parameter that we can tune. So using the same procedure as
before, we can iterate through 5 to 100 different clusters and see
what fits best.

K_range = range(knnmodel, :K, lower=5, upper=100);

knnTunedModel = TunedModel(model=knnmodel,
                           resampling = CV(nfolds=10, shuffle=true),
                           tuning = Grid(resolution=200),
                           range = K_range,
                           measures=[auc, log_loss, accuracy, kappa, brier_loss]);

knnTunedMachine = machine(knnTunedModel, X_trans, y);
fit!(knnTunedMachine, rows=train, verbosity=0)

Again, when we plot the results of the tuning we see the obvious
overfitting as the number of neighbours increases.

plot(knnTunedMachine)

report(knnTunedMachine).best_model

KNNClassifier(
  K = 9, 
  algorithm = :kdtree, 
  metric = Distances.Euclidean(0.0), 
  leafsize = 10, 
  reorder = true, 
  weights = NearestNeighborModels.Uniform())

9 clusters appear to be the optimal amount.

So no, a better log-loss and AUC measures, but the accuracy and
\(kappa\) metrics are worse.

Random Forest

Ok, final model! This is a random forest, so similar to the XGBoost
method.

randomforestModel = @load RandomForestClassifier pkg=DecisionTree verbosity = 0
randomforestmodel = randomforestModel()

rfMachine = machine(randomforestmodel, X_trans, y)
evaluate!(rfMachine,
     rows=train,
     resampling=CV(nfolds=6,shuffle=true),
     measures=[log_loss, accuracy, kappa, brier_loss, auc],
     verbosity=0)

It does very well with the best \(\kappa\) aside from the XGBoost
models.

Model Stacking

Ok, I lied, one more model, but this is a combination of all the
previous results. We have 4 different candidates and rather than
choosing the single best model, we can merge them to create a
super-model that blends the prediction. This is called model stacking
and comes in useful when the different model types perform well in
different conditions.

MLJ supports model stacking straight out the box, you just have to
pass it either the default or tuned model type. I pass the random
forest, elastic-net, knn, and XGBoost model into the stacking procedure
and see what comes out.

stackModel = Stack(;metalearner=logisticClassifier(lambda = 0),
                resampling=CV(nfolds = 6, shuffle=true),
                measures=[auc],
                rf = randomforestModel(),
                lm = report(lmTunedMachine).best_model,
                knn = report(knnTunedMachine).best_model,
                xgb = report(xgbTunedMachine).best_model)

stackedMachine = machine(stackModel, X_trans, y)

fit!(stackedMachine, rows=train, verbosity=0)

evaluate!(stackedMachine; 
          rows=train,
          resampling=CV(shuffle=true),  
          measures=[auc, accuracy, kappa, brier_loss, auc])

This stacked model does not perform as well as the XGBoost model on
its own, which is a shame.

That finishes all the model fitting. Given we only have 1000
observations, any more hyperparameter tuning is going to lead to
overfitting.

When we look at the accuracy and \(\kappa\) performance XGBoost comes
out the best and gives a 2.3% uplift on the null model. In the loan
context, this means being able to avoid an extra 2 defaulting loans
out of 100. The logistic regression models struggle to compare to the nonparametric methods, even with the penalisation in the elastic-net method.

The stacked model is slightly worse than the tuned XGBoost model. As we are just combining the different models this really highlights that the different model types aren’t capturing anything too different and the XGBoost model is sufficient on its own.

All of the above metrics are based on the training data though, so we will now evaluate the test data to see what model comes out on top.

Probability Calibration

We are relying on the probability of the prediction and whether it
reflects the true underlying probability of default on the loan. It’s
not enough to just predict whether the loan will default or not, we
have to get an idea of how good our probability output is. This is
where we assess how calibrated the model outputs are. Simply put, for
all the loans that we predict to have a 10% chance of defaulting, do they default at a rate of 10%? If they do, then we can say the
model is well-calibrated.

We will evaluate the calibration of our model on the test data and see how it lines up.

To produce a calibration plot we partition our predicted probabilities
into increasing groups and then calculate the number of loans that
each defaulted in those groups. So for the 0-10% bucket, selected all
the loans that we predicted a probability in that range and then
calculate how many defaulted. A good model will have around 10% of
them defaulting.

We do this for each model and plot the results.

xgboostProb = pdf.(MLJ.predict(xgbTunedMachine, rows=test), 1)
rfProb =pdf.(MLJ.predict(rfMachine, rows=test), 1)
knnProb = pdf.(MLJ.predict(knnTunedMachine, rows=test), 1)
stackProb = pdf.(MLJ.predict(stackedMachine, rows=test), 1)
lmProb = pdf.(MLJ.predict(lmTunedMachine, rows=test), 1)
irProb = pdf.(MLJ.predict(irMachine, rows=test), 1)

probFrame = DataFrame(BadLoan = Array(y[test]), 
                      XGBoost = xgboostProb, RandomForest = rfProb, 
                      KNN=knnProb, Stacked=stackProb,
					  LM = lmProb, IR = irProb)

probFrame = stack(probFrame, 2:7)

rename!(probFrame, ["BadLoan", "Model", "Prob"])

#Cut the probabilities into 10% groups.
lData = @transform(probFrame, :prob = cut(:Prob, (0:0.1:1.1)))
gData = groupby(lData, ["Model", "prob"])
calibData = @combine(gData, :N = length(:BadLoan), 
                            :DefaultRate = mean(:BadLoan), 
                            :PredictedProb = mean(:Prob))

calibData = @transform(calibData, :Err = 1.96 .* sqrt.((:PredictedProb .* (1 .- :PredictedProb)) ./ :N))


function calib_plot(calibData, m)

    p = plot(calibData[calibData.Model .== m, :PredictedProb], 
             calibData[calibData.Model .== m, :DefaultRate], 
             yerr = calibData[calibData.Model .== m, :Err],
        seriestype=:scatter, title=m, legend=:none, ylim=[0,1])
    plot!(p, 0.0:0.1:1.0, 0.0:0.1:1.0, label=:none)
end

plot(
    calib_plot(calibData, "KNN"), 
    calib_plot(calibData, "LM"),
    calib_plot(calibData, "XGBoost"),
    calib_plot(calibData, "Stacked"),
    calib_plot(calibData, "IR"),
    calib_plot(calibData, "RandomForest")
)

The stacked model goes a bit haywire, but overall the models and their confidence intervals line up with the perfectly calibrated red line.

We can now move on to the other metrics but this time apply them to
the test data.

modelNames = ["Null", "IR Only", "RF", "LM", "LM Tuned", "XGBoost", "XGBoost Tuned", "KNN", "KNN Tuned", "Stacked"]
modelMachines = [constMachine, irMachine, rfMachine,
               lmMachine, lmTunedMachine,
               xgbMachine, xgbTunedMachine,
               knnMachine, knnTunedMachine,
               stackedMachine]

With all the models in a list, we can map through the evaluation
metrics on the test data and see what the different models
produce. Any model that performed well on the training data but poorly
on the test data is overfitting.

aucRes = DataFrame(AUC = map(x->auc(MLJ.predict(x,rows=test), y[test]), 
                modelMachines),
          Model = modelNames)
kappaRes = DataFrame(Kappa = map(x->kappa(MLJ.predict_mode(x,rows=test), y[test]), modelMachines),
          Model = modelNames)
evalRes = leftjoin(aucRes, kappaRes, on =:Model)

So on the training dataset, we can see that the default XGBoost model
performs the best, indicating that tuning the hyperparameters just
leads to overfitting. As the dataset is so small this makes sense,
there needs to be an increase in the number of loans before we can
start seeing a performance benefit in adjusting the
hyperparameters.

To get a sense of all the model performances we can use a quadrant
plot of these two models to highlight how good they are.

evalResSub = @subset(evalRes, :Kappa .> 0)

plot(evalResSub.AUC, evalResSub.Kappa, seriestype = :scatter, group=evalResSub.Model, 
     legend=:none, series_annotations = text.(evalResSub.Model, :bottom, pointsize=8),
     xlabel = "AUC", ylabel = "Kappa")

So going forward, we will use the default XGBoost model as our
probability generator.

The confusion matrix helps some up the end goal of this project. We
want to invest in good loans and avoid the bad ones. Every good
loan that we don’t invest in because our model thought it would
default is an opportunity to make money lost, and likewise, every
loan that our model thought would be good that goes bad will cost us
money.

ConfusionMatrix()(MLJ.predict_mode(xgbMachine,rows=test), y[test])

              ┌───────────────────────────┐
              │       Ground Truth        │
┌─────────────┼─────────────┬─────────────┤
│  Predicted  │      0      │      1      │
├─────────────┼─────────────┼─────────────┤
│      0      │     626     │     29      │
├─────────────┼─────────────┼─────────────┤
│      1      │      4      │     39      │
└─────────────┴─────────────┴─────────────┘

• Opportunity lost -> predicted 1 but the truth was 0.
• Money lost -> predicted 0 but the truth was 1.

So in our test set we can that the lost opportunity is quite small,
the real danger is in the money lost category, all the times our model
predicted a loan wouldn’t default but it does.

Investing in Property Loans with Our Model

Right, we’ve got a model that we think is producing sensible
results. Now the important question is whether it helps us make money.

Looking at the test set we will predict whether it will
default and use that to inform our investment decision. If the loan
doesn’t default, we earn the interest rate, if it defaults, we lose
our invested capital. Side note, this is very much like sports betting
with asymmetrical payoffs.

Let’s only invest if there is less than a 50% chance that the loan defaults.

ir = X_encoded[test, "Interest_Rate"] ./100

investFrame = DataFrame(BadLoan = Array(y[test]), 
                        IR = ir, 
                        PBadLoan = pdf.(MLJ.predict(xgbMachine, rows=test), 1))

investFrame = @transform(investFrame, :Invest = Int.(:PBadLoan .<= 0.5))

unitInvest = @combine(groupby(investFrame, ["Invest", "BadLoan"]), 
         :N = length(:BadLoan), 
         :TheoReturn = sum(  (-1 .* (:BadLoan .== 1)) + (:BadLoan .== 0) .* :IR) )

So there are 43 loans we chose to not invest in (4+39). 4 of those went on to pay back their loan, so an opportunity cost of 0.43. But 39 did default, so we avoided a big miss there.
We invested in 626 loans, of which 29 defaulted. But our profits on
the good loans outweighed these losses. So overall, happy days, we
ended up with more money than we started.

Our 50% threshold for investing was arbitrary, so we can vary it and see how the profit changes.

function thresh_experiment(investFrame, thresh)
    investFrame = @transform(investFrame, :Invest = Int.(:PBadLoan .<= thresh), :Thresh = thresh)

    @combine(groupby(investFrame, ["Thresh", "Invest", "BadLoan"]), 
         :N = length(:BadLoan), 
         :TheoReturn = sum(  (-1 .* (:BadLoan .== 1)) + (:BadLoan .== 0) .* :IR) )
end
threshRes = vcat(thresh_experiment.([investFrame], 0.05:0.05:0.95)...);

profitRes = @combine(
    groupby(
    @subset(threshRes, :Invest .== 1),
    "Thresh"
    ),
    :Profit = sum(:TheoReturn)
)

plot(profitRes.Thresh, profitRes.Profit, xlabel="Threshold Probability", ylabel="Profit", label = :none)

We get to improve the profitability slightly, but not in a meaningful
way, so let’s stick to our 50% rule. Why complicate things?

Kelly Betting

The Kelly bet is the optimal position sizing based on the expected
payoff and independent estimation of the probability that payoff
happens. This is perfectly suited to our situation and allows us to
invest more into loans where there is a dislocation between its
interest rate and probability of default.

From wikipedia the Kelly bet formula is:

\[f = p + \frac{p-1}{b},\]

where \(p\) is the probability the bet comes off, \(b\) is the proportion of the bet won.

In our case, the bet is place 1 unit and return \(1+r\) units. So \(b = (1+r-1)/1\) which comes our to \(r\).

So the formula is:

\[\begin{align}
f & = (1-p_\text{default}) + (1-p_\text{default} – 1)/r, \\
& = 1 – p_\text{default} – p_\text{default}/r, \\
& = 1 – p_\text{default}(1 – 1/r),
\end{align}\]

we can calculate this bet size for each loan and see how it changes our profitability.

investFrame = @transform(investFrame, :KellyBet = max.((1 .- :PBadLoan) .- (:PBadLoan./:IR), 0))
kellyInvest = @combine(groupby(investFrame, ["Invest", "BadLoan"]), 
         :N = length(:BadLoan), 
         :TotalStaked = sum(:KellyBet),
         :Profit = sum(  (-:KellyBet .* (:BadLoan .== 1)) + (:BadLoan .== 0) .* :KellyBet .* :IR) )

kellyVsUnit = @combine(investFrame, 
         :N = length(:BadLoan), 
         :TotalKellyStaked = sum(:KellyBet),
         :KellyProfit = sum(  (-:KellyBet .* (:BadLoan .== 1)) + (:BadLoan .== 0) .* :KellyBet .* :IR),
         :TotalUnitStaked = sum(:Invest),
         :UnitProfit = sum( :Invest .* ( (-1 .* (:BadLoan .== 1)) + (:BadLoan .== 0) .* :IR)))
@transform(kellyVsUnit, :KellyROIC = :KellyProfit ./ :TotalKellyStaked, 
                        :UnitROIC = :UnitProfit ./ :TotalUnitStaked)

This shows that Kelly betting helps manage the risk of the
strategy and accounts for the uncertainty in the prediction. It has
about a 1% improvement in return versus just betting one unit each
time. It’s not the true return, as it isn’t reinvesting the profits from the
earlier paid-back loans, but it is a good indication of the
profitability of the system.

Where Will it Go Wrong?

Like all investment strategies we have to think about how we might end
up with our faces ripped off. The first one is that I cannot
understate the risk of investing in bridge loans for properties in
Europe. This is highly speculative and there is probably a reason that
the companies are coming to retail investors rather than other sources
of finance. So, there will be an element of selection bias in these
loans, they will perform well until they don’t.

Next up is the macro environment. Interest rates are currently rising
and the world is moving into a new regime. We can’t imagine that our
model will account for this so the underlying default is likely to
change over the next few months. As credit conditions tighten, it’s
likely these loans and new upcoming loans have a higher rate of
default.

Estateguru could just up and leave too. There is significant
counterparty risk in this trade as Estateguru is hardly an
established company. They recently raised money on Seedrs, so still a
very early-stage company.

Lack of sample size. This model has only 3789 rows of data of which 2327 can be used to learn about the default probability. So in reality we are very uncertain.

Currency risk. All the loans are denominated in EUR and as I get paid in GBP I’ll have to hedge out the fluctuations in currency. Can’t have the high-risk returns being eaten away by the UKs monetary policy.

Plus the most dangerous unknowns, the ones we haven’t got a clue about. We didn’t see COVID coming, who knows what else could be on the horizon.

Still, we can at least look at the current open loans and see what the model would say.

Open Loans

These are loans that we can invest in today if so inclined. So let’s estimate their probability of default.

I’ll apply all the previous transformations and predict the model

openLoans[!, "Property_Type_A"] .= "N/A"
openLoans[!, "Property_Type_B"] .= "N/A"

for i in 1:nrow(openLoans)

    pt = strip.(split(openLoans[i, "Property_Type"], "-"))
    openLoans[i, "Property_Type_A"] = pt[1]
    
    if length(pt) == 1
        openLoans[i, "Property_Type_B"] = pt[1]
    else
        openLoans[i, "Property_Type_B"] = join(pt[2:end], " ")
    end
end

@combine(groupby(openLoans, ["Property_Type_A", "Property_Type_B"]), :n=length(:Status))

Not very interesting property types! At least they are all the same
category as the training set.

openData = coerce(openLoans,
                      :Country=>Multiclass,
                      :Schedule_Type => Multiclass,
                      :Property_Type_A => Multiclass,
                      :Property_Type_B => Multiclass,
                      :Loan_Type => Multiclass,
                      :Suretyship_existence => Multiclass);


openData[!, "log_Funded_Total_Amount"] = log.(openData[!, "Funded_Total_Amount"])
openData[!, "log_Property_Value"] = log.(openData[!, "Property_Value"])
openData[!, "LTV"] = openData[!, "LTV"] /100
openData[!, "Loan_Period"] = openData[!, "Loan_Period"] /12

openData = select(openData, Not([:Funded_Total_Amount, :Property_Value]));

for col in ["Country", "Schedule_Type", "Property_Type_A", "Property_Type_B", 
            "Loan_Type", "Suretyship_existence"]
    levels!(openData[!, col], levels(modelData2[!, col]))
end


X_encoded_open = MLJ.transform(encMach, openData)
X_trans_open = MLJ.transform(stanMach, X_encoded_open)

pDefault = pdf.(MLJ.predict(xgbMachine, X_trans_open), 1)
openPred = DataFrame(LoanID = openLoans.Loan_code, IR = openLoans[!, "Interest_Rate"]./100, 
                     PBadLoan = pDefault)
@transform(openPred, :KellyBet = max.((1 .- :PBadLoan) .- (:PBadLoan./:IR), 0))

So six properties should be avoided, (the ones where the Kelly bet is
0) and the model believes all others can be invested in. So an encouraging result.

Conclusion

Given the open data, we have built a nice model exploring the
probability of default and turned it into a profitable investing
strategy. You hopefully know a bit more about MLJ.jl and how it can be
your one-stop-shop for machine learning in Julia.

Now like most things in this alt-finance space it is incredibly risky
to invest in these loans and given the way interest rates are going
it’s likely that the default rate will skyrocket. So don’t go
investing in these loans because of what I said, I’m no expert.

By: Dean Markwick's Blog -- Julia

Re-posted from: https://dm13450.github.io/2022/05/11/modelling-microstructure-noise-using-hawkes-processes.html

Microstructure noise is where the
price we observe isn’t the ‘true’ price of the underlying asset. The
observed price doesn’t diffuse as we assume in your typical
derivative pricing models, but instead, we see
some quirks in the underlying data. For example, there is an explosion of
realised variance as we use finer and finer time subsampling periods.

Last month I wrote about
calculating realised volatility
and now I’ll be taking it a step further. I’ll show you how this
microstructure noise manifests itself in futures trade data and how I
use a Hawkes process to come up with a price formation model that fits
the actual data.

The original work was all done in
Modeling microstructure noise with mutually exciting point processes. I’ll
be explaining the maths behind it, showing you how to fit the models
in Julia and hopefully educating you on this high-frequency finance
topic.

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!

A Bit of Background

My Ph.D. was all about using Hawkes processes to model when things
happen and how these things are self-exciting. You may have read my
work on Hawkes process before, either through my Julia package
HawkesProcess.jl,
or examples of me
using Hawkes processes to model terror attack data or just
how to calculate the deviance information criteria of a Hawkes process. The
real hardcore might have read my Ph.D. thesis.

But how can we use a Hawkes process to model the price of some
financial asset? There is a vast amount of research and work about how Hawkes
processes can be used in price formation processes for financial
instruments and high-frequency types of problems and this post will
act as a primer to anyone interested in using Hawkes processes for
such problems.

The High-Frequency Problem

At short timescales, a price moves randomly rather than with any
trend. Amazon stock might trade thousands of times in a minute but
that’s supply and demand changing rather than people thinking
Amazon’s future is going to change from one minute to the next. So we need
a different way of thinking about how prices move at short timescales
compared to longer timescales.

We can build nice mathematical models guessing how a price might move;
it might move like a random walk or maybe a random walk with jumps in
the price now and then. But, no matter what model we use, it must
match up with what is observed in the real world. One of
these observations is a phenomenon called ‘microstructure noise’ and we
only see it in high-frequency data.

Microstructure noise is a catch-all term for different things happening
in the market. This includes, bid-ask bounce, people buy and sell at
two different prices, so looks like the price is moving, but in
reality, just oscillating around the mid-price. The discreteness of
prices at these time scales also plays a part. There is a minimum
increment level that prices change buy and this can have real effects
on how prices move. Exchanges need to pay attention to their tick
sizes, as it can help or hinder liquidity if they are set
incorrectly. This then has a real effect that we can observe when we
calculate a realised volatility.

Realised volatility is measuring how much the price moved in a
period. These high-frequency effects are going to give the impression
of more movement than what the ‘true’ volatility is. So we end up
seeing our measurement of volatility explode as the time scale we use
gets smaller and smaller. Calculating the volatility using 1-second
intervals gives a larger value than if we used 1 minute
intervals. This means that our volatility estimation depends on the
time scale used, so what is the ‘real’ volatility?

We aren’t interested in working out what the real volatility
is. Instead, we want to build a model for a price that displays this
volatility scaling effect.

The Hawkes Process Model for a Price

How do we use Hawkes processes to build a model that will have this microstructure noise?

Lets call the price at time \(t\), \(X(t)\) and guess that it moves by
summing up the positive jumps \(N_1(t)\) and negative jumps \(N_2(t)\)
that also happen at time \(t\).

\[X(t) = N_1(t) – N_2(t)\]

When do these jumps occur? How are these jumps distributed? This is
where we use the Hawkes process.

A Hawkes process is a self-exciting point process. When something
happens it increases the probability of another event happening. This
is the self-exciting behaviour we want. Every time there is a positive
jump, there is an increase in the probability of a negative jump
happening and, likewise, when there is a negative jump there is a greater
probability of a positive jump.

Each jump causes the probability of the other jump happening like in
this picture

When someone buys and pushes the price higher by removing that
liquidity, it’s more likely that someone will now sell at the new
higher price introducing some downward pressure. This is mean
reversion where prices move higher and then outside forces
push it lower and vice versa.

With Hawkes processes there are three parameters:

• The background rate or when the jumps randomly occur. This is common to both positive and negative jumps.
• \(\kappa\) – the ‘force’ that pushes and increases the probability of the other jump happening.
• The kernel \(g(t)\) dictates how long the force lasts. This is an exponential decay with parameter \(\beta\).

We will fit a Hawkes process to a price series
to infer the 3 parameters that describe how the jumps behave. This
model will hopefully replicate the ‘microstructure noise’ effects we
see in practise.

Futures Trade Data

In the early days of my Ph.D., I answered an email that was advertising
for early grad students to do some prop trading. As part of the
interview, they gave me some data and asked me to write a simple
moving cross-over strategy. I failed miserably as I never
heard back from them. But I did get some nice data, which now that I’m
older and wiser, recognise as reported trades from a futures
exchange. This is the data we will be using today to calculate the
mode and luckily it’s similar to the data they use in the original
paper.

using CSV
using DataFrames, DataFramesMeta
using Plots
using Dates
using Statistics

All the usual packages when working with data in Julia.

rawData = CSV.read("fgbl__BNH14_clean.csv", DataFrame, header=false)
rename!(rawData, [:UnixTime, :Price, :Volume, :DateTime]);
first(rawData, 5)

To clean the data we convert the Unix timestamp to an actual DateTime object and pull out the hour and the date of the trade.

cleanData = @transform(rawData, DateTimeClean = DateTime.(:DateTime, dateformat"mm/dd/yyyyHH:MM:SS.sss"), 
                                DateTimeUnix = unix2datetime.(:UnixTime ./ 1000) )
cleanData = @transform(cleanData, Date = Date.(:DateTimeUnix),
                                  Hour = hour.(:DateTimeUnix));
first(cleanData[:,[:UnixTime, :DateTimeClean, :DateTimeUnix]], 5)

To get an idea of the data we are looking at I aggregate the total
number of trades and total volume of the trades over each day and plot
that as a time series.

dayData = groupby(cleanData, :Date)
dailyVolumes = @combine(dayData, TotalVolume = sum(:Volume),
                                  TotalTrades = length(:Volume),
                                  FirstTradeTime = minimum(:DateTimeUnix))

xticks = minimum(dailyVolumes.Date):Month(2):maximum(dailyVolumes.Date)
xticks_labels = Dates.format.(xticks, "yyyy-mm")

p1 = plot(dailyVolumes.Date, dailyVolumes.TotalVolume, seriestype=:scatter, label="Daily Volume", legend = :topleft, xticks = (xticks, xticks_labels))
p2 = plot(dailyVolumes.Date, dailyVolumes.TotalTrades, seriestype=:scatter, label= "Daily Number of Trades", legend = :topleft, xticks = (xticks, xticks_labels))
plot(p1, p2, fmt=:png)

It takes a while for the trading to take off in this future
contract. This is where it slowly becomes the front-month contract and
then is the most active.

Also, because trading doesn’t cross over the daylight saving dates, we don’t have to worry about timezones. Always a bonus!

What about if we look at what hour is the most active?

hourDataG = groupby(cleanData, [:Date, :Hour])
hourDataS = @combine(hourDataG, TotalHourVolume = sum(:Volume),
                                  TotalHourTrades = length(:Volume))
hourDataS = leftjoin(hourDataS, dailyVolumes, on=:Date)

hourDataS = @transform(hourDataS, FracVolume = :TotalHourVolume ./ :TotalVolume)

hourDataG = groupby(hourDataS, :Hour)
hourDataS = @combine(hourDataG, MeanFracVolume = mean(:FracVolume),
                                 MedianFracVolume = median(:FracVolume))
sort!(hourDataS, :Hour)
bar(hourDataS.Hour, hourDataS.MedianFracVolume * 100, title = "Fraction of Daily Volume", label=:none, fmt=:png)

Early in the morning (just after the exchange opens) and late afternoon (when the Americans start trading) is where there is the most activity.

For our analysis, we are going to be focused on the hours 14, 15, 16
to make sure that we have the most active period and this is the same
as what the original paper did, took a subset of the day.

How to Calculate the Volatility Signature

Let \(X(t)\) be the price of the future at time \(t\). The signature is the quadratic variation over a window of \([0, T]\), which is more commonly known as the realised volatility:

\[C(\tau) = \frac{1}{T} \Sigma _{n=0} ^{T/\tau} \mid X((n+1) \tau) – X(n \tau) \mid ^2 .\]

\(\tau\) is our sampling frequency, say every minute, etc.

To calculate the volatility across the trades we have to pay particular attention to the fact that these trades are irregularly spaced, so we need to fill forward the price for every \(t\) value.

function get_price(t::Number, prices, times)
    ind = min(searchsortedfirst(times, t), length(times))
    sp = ind == 0 ? 0 : prices[ind]
end

function get_price(t::Array{<:Number}, prices, times)
    res = Array{Float64}(undef, length(t))
    for i in eachindex(t)
        res[i] = get_price(t[i], prices, times)
    end
    res
end

Our get_price function here will return the last price before time \(t\).

To calculate the signature value we chose a \(\tau\) value, generate
the indexes between 0 and \(T\) using a \(\tau\) step size. Pull the
price at those times and calculate the quadratic variation. Again, we
add a method to calculate the signature for different \(\tau\)’s.

function signature(tau::Number, x, t, maxT)
    inds = collect(0:tau:maxT)
    prices = get_price(inds, x, t)
    
    rets = prices[2:end] .- prices[1:(end-1)]
    (1/maxT) * sum(abs.(rets) .^ 2)
end

function signature(tau::Array{<:Number}, x, t, maxT)
   res = Array{Float64}(undef, length(tau))
    for i in eachindex(res)
        res[i] = signature(tau[i], x, t, maxT)
    end
    res
end

Now let’s apply this to the data. We are only interested when the
future was actively trading and in the hours between 14:00 and 16:00.
We convert the times into seconds since 15:59:59 and calculate the
signature for all the dates, before taking the final average.

We are taking the log price of the last trade to represent our actual
\(X(t)\). We just look at 2014 dates too as that is when the future is
active.

cleanData2014 = @subset(cleanData, :Date .>= Date("2014-01-01"))

uniqueDates = unique(cleanData2014.Date)

eventList = Array{Array{Float64}}(undef, length(uniqueDates))
priceList = Array{Array{Float64}}(undef, length(uniqueDates))

signatureList = Array{Array{Float64}}(undef, length(uniqueDates))
avgSignature = zeros(length(1:1:200))

for (i, dt) in enumerate(uniqueDates)
   
    subData = @subset(cleanData2014, :Date .== dt, :Hour .<= 16, :Hour .>= 14)
    eventList[i] = getfield.(subData.DateTimeClean .- (DateTime(dt) + Hour(14) - Second(1)), :value) ./ 1e3
    priceList[i] = subData.Price
    
    signatureList[i] = signature(collect(1:1:200), log.(priceList[i]), eventList[i] .+ rand(length(eventList[i])), 3*60*60 + 1)
    avgSignature .+= signatureList[i]
end

avgSignature = avgSignature ./ length(eventList);

To plot the signature we take the average across all the dates and
then normalised by the \(\tau = 60\) value.

plot(avgSignature / avgSignature[60], seriestype=:scatter, 
    label = "Average Signature", xlabel = "Tau", ylabel = "Realised Volatility (normalised)", fmt=:png)

This is an interesting plot with big consequences in high-frequency finance.

This explosion in realised volatility at small timescales (\(\tau
\rightarrow 0\)) comes from microstructure noise. If prices evolved
as Brownian motion, the above plot would be flat for all timescales so
the above result contradicts lots of classical finance assumptions.

Practically, this is a pain if we are trying to measure the currently
volatility, it depends on the timescale we are looking at, there isn’t
one true volatility using the normal methods. Instead, we need to be
aware of these microstructure effects as we use a finer \(\tau\) over
which to calculate the volatility.

This is where the Hawkes model comes in. If we assume the price,
\(X(t)\) moves as stated in Equation (), can we produce a similar signature plot?

The Theoretical Signature Under a Hawkes Process

After doing some maths (you can read the paper for the full details), we arrive at the following equation for the theoretical signature.

If both \(N_1\) and \(N_2\) are realisations of Hawkes processes with
parameters \(\mu, \kappa\) and \(g(t) = \beta \exp (-\beta t)\) then their intensity can be written as

\[C(\tau) = \Lambda \left( k ^2 + (1 – k ^2) \frac{1 – e ^{-\gamma \tau}}{ \gamma \tau} \right),\]

where

\(\Lambda = \frac{2 \mu}{1 – \kappa}, k = \frac{1}{1 + \kappa}, \gamma = \beta (\kappa + 1)\).

These are from the paper and adjusted based on my parameterisation of the Hawkes process. This gives us our theo_signature function.

function theo_signature(tau, bg, kappa, kern)
    Lambda = 2*bg/(1-kappa)
    k = 1/(1 + kappa)
    gamma = kern*(kappa + 1)
    @. Lambda * (k^2 + (1-k^2) * (1 - exp(-gamma * tau)) / (gamma*tau))
end

Calibrating the Hawkes Process Model

We now move on to fitting the Hawkes process to the data. I’ll be using
a new method that takes a different approach to my package HawkesProcesses.jl.

We have a theoretical volatility signature from a Hawkes process
(theo_signature) and the above plot of what the actual signature
looks like it. Therefore, it is just a case of optimising over a loss
function to find the best fitting parameters. I’ll use root mean
square error as my loss function and simply use the Optim.jl package
to perform the minimisation.

signatureRMSE(x, sig) = sqrt(
    mean(
        (sig .- theo_signature(1:200, x[1], x[2], x[3])).^2
        )
)

using Optim
optRes = optimize(x->signatureRMSE(x, avgSignature/avgSignature[10]), rand(3))

paramEst = Optim.minimizer(optRes)

3-element Vector{Float64}:
 0.24402236592012655
 0.7417867072115396
 0.19569169443936185

These are the three parameters of the Hawkes process which appear
sensible.

plot(avgSignature, label="Observed", seriestype=:scatter)
plot!(avgSignature[10]*theo_signature(1:200, paramEst[1], paramEst[2], paramEst[3]), label="Theoretical", lw=3, xlabel = "Tau", ylabel = "Realised Variance", fmt=:png)

So a nice match-up between the theoretical signature and what we
observed. This gives some weight to the Hawkes model as a
representation of the price process.

Interpreting the Hawkes Parameters

Our \(\kappa\) value of 0.75 shows there is a large amount of
excitement with each price jump. a \(\beta\) value of 0.2 shows that this
mean reversion lasts for about 5 seconds. So if we see an uptick in
the price, we expect a downtick with a rough half-life of five
seconds. The opposite is also true, a downtick likely leads to an
uptick 5 seconds later.

Conclusion

Microstructure noise shows up when we start calculating volatility on
a high-frequency timescale. We have shown that it is a real effect
using some futures data and then built a Hawkes model to try and
reproduce this effect. We managed to get the right shape of the
volatility signature and found that was quite a bit of mean reversion
between the up and downticks that lasted around 5 seconds.

What’s next? In a future post, I will introduce another dimension and
show and there can also be correlation across assets under a similar
method to reproduce another high-frequency phenomenon. This will be
based on the same paper and show you how we can start looking at the
correlation between two assets and how that changes at high-frequency
time scales.