U NIVERSITY OF C APE T OWN

H ONOUR ’ S P ROJECT R EPORT

An Exercise in R: High Frequency Covariance estimation using

Malliavin-Mancino and Hayashi-Yoshida estimators

Author: Supervisor:
Patrick C HANG & Roger Assoc. Prof. Tim G EBBIE
B UKURU

A report submitted in fulfilment of the requirements

for the degree of Honours in Statistical Science
in the

Department of Statistical Sciences

October 21, 2019

 iii

Declaration of Authorship
We, Patrick C HANG & Roger B UKURU , declare that this report titled, “An Ex-
ercise in R: High Frequency Covariance estimation using Malliavin-Mancino
and Hayashi-Yoshida estimators” and the work presented in it are our own.
We confirm that:
• This work was done wholly or mainly while in candidature for a re-
search degree at this University.
• Where any part of this thesis has previously been submitted for a de-
gree or any other qualification at this University or any other institu-
tion, this has been clearly stated.
• Where we have consulted the published work of others, this is always
clearly attributed.
• Where we have quoted from the work of others, the source is always
given. With the exception of such quotations, this thesis is entirely our
own work.
• We have acknowledged all main sources of help.
• Where the thesis is based on work done by ourselves jointly with oth-
ers, we have made clear exactly what was done by others and what we
have contributed ourselves.
Signed:

Date:
 v

UNIVERSITY OF CAPE TOWN

Abstract
Science
Department of Statistical Sciences

Honours in Statistical Science

An Exercise in R: High Frequency Covariance estimation using

Malliavin-Mancino and Hayashi-Yoshida estimators

by Patrick C HANG & Roger B UKURU

We revisit and demonstrate two well-known non-parametric estimators; the

Malliavin-Mancino (MM) [1], [2] and the Hayashi-Yoshida (HY) [3] estima-
tor. Both address the issue of covariance estimation for high-frequency asyn-
chronous time-series data. The first by embracing a Fourier perspective and
the latter using averaging over discrete windows. The aim of the work here
aims to provide an easily re-used practical tutorial to argue for the efficacy of
the MM estimator in high-frequency finance applications. Towards this end,
we conduct Monte Carlo experiments to demonstrate that the two estima-
tors differ only under asynchronous observations, where the MM estimator
has lower correlation estimates compared to the HY estimator. Unsurpris-
ingly, we attribute this difference to the Epps effect [4]. However, as a novel
application, we show the existence of the Epps effect in the top 10 stocks
from the Johannesburg Stock Exchange (JSE) by various methods of aggre-
gating Trade and Quote (TAQ) data. Specifically, by comparing calendar
time based sampling with volume time sampling methods. We argue that
the MM estimator is more representative of trade-time reality. The world
of high-frequency finance is not a missing data problem as sampled from
some underlying synchronous continuous stochastic process, but rather it is a
world of truly disconnected, inter-related, discrete and asynchronous events
where the relationship between events are the fundamental entities and mea-
surables. We argue the MM estimator is a more faithful representation of the
underlying data as it does not over-estimate short-term correlations in such
an asynchronous event driven world.
 vii

Acknowledgements
We would firstly like to thank our supervisor Tim Gebbie for guiding us
through the various difficulties of this project and never revealing the answer
straight away, giving us the satisfaction of solving the problem ourselves. Ef-
fectively guiding us to ensure we were on time for the deadlines and always
finding the time to answer our questions.
We would secondly like to thank Etienne Pienaar and Melusi Mavuso for
their input and assistance. We would also like to thank Diane Wilcox, Chanel
Malherbe and Dieter Hendricks who provided us with resources to fast track
our progress.
Patrick Chang would like to thank the National Research Foundation for
funding his study at the University of Cape Town.
Thank You
Patrick C HANG & Roger B UKURU
 ix

Contents

Declaration of Authorship iii

Abstract v

Acknowledgements vii

1 Introduction 1

2 The Theory 3
2.1 Trigonometric Fourier Transform . . . . . . . . . . . . . . . . . 3
2.1.1 Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1.2 Numerical Implementation . . . . . . . . . . . . . . . . 5
2.2 Complex Exponential Fourier Transform . . . . . . . . . . . . . 6
2.2.1 Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2.2 Numerical Implementation . . . . . . . . . . . . . . . . 8
2.3 Hayashi-Yoshida Estimator . . . . . . . . . . . . . . . . . . . . 10
2.3.1 Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3.2 Numerical Implementation . . . . . . . . . . . . . . . . 11

3 Monte Carlo Experiments 13

3.1 Effect of missing observations . . . . . . . . . . . . . . . . . . . 13
3.2 Effect of the SDE . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3 Effect of Asynchrony . . . . . . . . . . . . . . . . . . . . . . . . 19

4 Data Engineering 25
4.1 Data Collection . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.1.1 Collecting Data From Bloomberg . . . . . . . . . . . . . 25
4.1.2 Updated HF-Data Pipeline . . . . . . . . . . . . . . . . 26
4.2 Data Cleaning . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.2.1 Data Types . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.2.2 Aggregation . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.2.3 Overnight returns . . . . . . . . . . . . . . . . . . . . . . 29
4.3 Creating Data Samples . . . . . . . . . . . . . . . . . . . . . . . 30
4.3.1 Asynchronous Data . . . . . . . . . . . . . . . . . . . . 30
4.3.2 Calendar Time TAQ Data Aggregation . . . . . . . . . 31
4.3.3 Intrinsic Time TAQ Data Aggregation . . . . . . . . . . 34
Derman Framework . . . . . . . . . . . . . . . . . . . . 34
Lining Up Events . . . . . . . . . . . . . . . . . . . . . . 36

5 Data Science 39
x

5.1 Calendar Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

5.1.1 Closing Prices . . . . . . . . . . . . . . . . . . . . . . . . 39
5.1.2 Volume Weighted Prices . . . . . . . . . . . . . . . . . . 41
5.2 Intrinsic Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.2.1 Derman Framework . . . . . . . . . . . . . . . . . . . . 44
5.2.2 Lining Up Events . . . . . . . . . . . . . . . . . . . . . . 47

6 Concluding Remarks 49

A Supporting Algorithms 51

B Appendix Derivation 55
B.1 Proof for Theorem 2.1.1 and 2.1.2 . . . . . . . . . . . . . . . . . 55
B.1.1 Proof for aq (Σ) . . . . . . . . . . . . . . . . . . . . . . . 55
B.1.2 Proof for a0 (Σ) . . . . . . . . . . . . . . . . . . . . . . . 59
B.1.3 Proof for bq (Σ) . . . . . . . . . . . . . . . . . . . . . . . 60
B.1.4 Proof for aq (Σi,j ) . . . . . . . . . . . . . . . . . . . . . . 61
B.1.5 Proof for a0 (Σi,j ) . . . . . . . . . . . . . . . . . . . . . . 62
B.1.6 Proof for bq (Σi,j ) . . . . . . . . . . . . . . . . . . . . . . 62
B.2 Proof for Theorem 2.2.1 . . . . . . . . . . . . . . . . . . . . . . . 63
B.3 Proof for Theorem 2.3.1 . . . . . . . . . . . . . . . . . . . . . . . 65

C Miscellaneous 73
C.1 Run times and General issues . . . . . . . . . . . . . . . . . . . 73
C.2 Supporting plots . . . . . . . . . . . . . . . . . . . . . . . . . . 75
 xi

List of Figures

3.1 The effect of missing data. . . . . . . . . . . . . . . . . . . . . . 14

3.2 The effect of various diffusion processes. . . . . . . . . . . . . . 16
3.3 The effect of volatility clustering and mean-reversion. . . . . . 18
3.4 Recovering the results from [12]. . . . . . . . . . . . . . . . . . 20
3.5 The effect of asynchrony. . . . . . . . . . . . . . . . . . . . . . . 22

4.1 Data Collection Pipeline . . . . . . . . . . . . . . . . . . . . . . 26

4.2 Trades with the same time stamp. . . . . . . . . . . . . . . . . . 28
4.3 Trades aggregated into a unique time stamp. . . . . . . . . . . 29
4.4 BTI and NPN Asynchronous Price Sample . . . . . . . . . . . 30
4.5 BTI and NPN Asynchronous Return Sample . . . . . . . . . . 31
4.6 BTI and NPN 10 Minute Closing Bar Return Sample . . . . . . 33
4.7 BTI and NPN 10 Minute VWAP Bar Return Sample . . . . . . 33
4.8 Multiple Ticker Derman Volume Buckets . . . . . . . . . . . . 35
4.9 Multiple Ticker Lining Up Events Volume Buckets . . . . . . . 37

5.1 Comparing the two estimators with closing price aggregation. 40

5.2 Comparing the two estimators with VWAP aggregation. . . . 43
5.3 Comparing the two estimators using Derman’s [14] intrinsic
time sampling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.4 Comparing the two estimators using a sampling method to
line up events. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

C.1 Comparing run times of various algorithms. . . . . . . . . . . 73

C.2 Price paths, returns and QQ-plot for the various SDEs. . . . . 75
C.3 Comparing the two estimators with the mistaken closing price
aggregation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
C.4 Comparing the two estimators with the mistaken VWAP ag-
gregation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
 xiii

List of Tables

4.1 Bloomberg APIs compared . . . . . . . . . . . . . . . . . . . . . 25

 xv

List of Abbreviations

API Application Programming Interface

FFT Fast Fourier Transform
GARCH Generalized Autoregressive Conditional Heteroskedasticity
GBM Geometric Brownian Motion
HY Hayashi Yoshida
JSE Johannesburg Stock Exchange
MM Malliavin Mancino
NaN Not a Number
NUFFT Non Uniform Fast Fourier Transform
OHCLV Open High Low Close Volume
OU Ornstein Uhlenbeck
RV Realised Volatility
SDE Stochastic Differential Equation
TAQ Trade And Qoute
VG Variance Gamma
VWAP Volume Weighted Average Price
 1

Chapter 1

Introduction

Covariation is a key parameter in finance with traditional applications in

portfolio optimisation and more recent applications in unsupervised state
discovery to discern changes in the system [5]. The availability of high-
frequency financial data has allowed the miss-estimation of large portfolio
correlation measures to be elevated by removing the problem of scarcity
of data [6]. However, having high-frequency financial data comes with its
caveats, specifically estimating the correlation becomes a much harder task
due to the asynchrony arising from having tick-by-tick trade data. Thus the
popular approach using the realised volatility estimator is problematic under
high-frequency asynchronous data, as it requires a choice of synchronisation
and data interpolation which leads to biases induced in the estimate [3], [7].
In this report, we present two non-parametric estimators designed specifi-
cally to deal with the asynchrony. The first estimator proposed by Malliavin
and Mancino [1], [2], [8] adopts a Fourier approach and the second estima-
tor proposed by Hayashi and Yoshida [3] uses the contributions from over-
lapping intervals to overcome the problem faced by the traditional realised
volatility.
We will investigate the correlation under high-frequency finance using the
aforementioned estimators in a data-informed approach. This novel approach
is due to how the price process is generated in the financial market - we ar-
gue that the price generation is better represented as discrete, asynchronous
events rather than samples from an underlying continuous stochastic pro-
cess; thus we abstain from taking the more popular market microstructure
noise framework [9]–[11]. A natural consequence of studying high-frequency
financial data is the need to address the Epps effect [4] which is the drop in
correlation associated with smaller sampling intervals. We demonstrate the
Epps effect arising from asynchrony, specifically the relation this has with the
level of asynchrony and the sampling interval under consideration [6], [12],
[13], furthermore we will investigate the Epps effect from a new perspective,
specifically by aggregating Trade and Quote (TAQ) data with different tech-
niques.
To this end, the report is structured as follows: in Section 2 we will cover the
2 Chapter 1. Introduction

brief derivation and implementation algorithms for the MM and HY estima-

tors. Section 3 we conduct Monte Carlo experiments in an attempt to iden-
tify how, when and where the two estimators differ. Therefore we begin the
comparison by assessing how the two estimators perform under asynchrony
- specifically asynchrony induced from a missing data manner. We then com-
pare the two estimators with various stochastic processes. Finally, we recover
the results from [12] and adapt his experiment to highlight the differences in
our data-informed approach compared to the market microstructure noise
approach. Section 4 outlines the data collection process and the various algo-
rithms employed to aggregate TAQ data from different perspectives, specifi-
cally we consider Closing and Volume Weighted Average Price (VWAP) bars
from the Calendar time approach. In addition, we consider an Intrinsic time
approach for which we will employ the framework provided by Derman [14]
and our method of aggregation in Intrinsic time. Section 5 combines the vari-
ous aggregation techniques from Section 4 along with the estimators to study
the Epps effect in the Johannesburg Stock Exchange (JSE). Finally, we end off
with Section 6 where we discuss some future topics of investigation and we
highlight the aspects of what we could and could not achieve with this study.
 3

Chapter 2

The Theory

2.1 Trigonometric Fourier Transform

2.1.1 Derivation
We present the derivation proposed by Malliavin and Mancino [1] which is
an estimator that is constructed in the frequency domain. The only assump-
tion required for the derivation is the Bachelier hypothesis. Which states that
all measurable economic data p∗ are driven by semi martingales which can
be decomposed into a drift term with bounded variation paths and a local
martingale [1], [15]. Thus their Itô stochastic differential equation given by

d
∑ σi dW i + β j dt.
j
dp j = (2.1.1)
i =1

If we denote Si (t) to be the generic asset price at time t, we will set pi (t) =
ln(Si (t)). It can be shown that the covariance matrix of diffusion processes
given by (2.1.1) can be presented as

d
∑ σi (t)σik (t).
j,k j
Σ (t) = (2.1.2)
i =1

For simplification purposes, we can always reduce a semi-martingale on a

fixed time window to the case where the window is [0, 2π ] by the change of
origin and re-scaling the unit of time [1]. We now define the Fourier coeffi-
cients of dp j as

1 2π Z
j
a0 (dp ) = dp j (t),
2π 0
1 2π
Z
j
ak (dp ) = cos(kt)dp j (t), (2.1.3)
π 0
1 2π
Z
j
bk (dp ) = sin(kt)dp j (t),
π 0

∀k ≥ 1. Similarly, the Fourier coefficients of the volatility is defined as

4 Chapter 2. The Theory

1 2πZ
a0 ( Σ ) = Σ(t)dt,
2π 0
1 2π
Z
ak (Σ) = cos(kt)Σ(t)dt, (2.1.4)
π 0
1 2π
Z
bk ( Σ ) = sin(kt)Σ(t)dt.
π 0

The main idea behind this method is to find a mathematical expression of the
Fourier coefficients of Σ using the Fourier coefficients of dp j [1]. This leads to
Theorem 3.1 in [1]:
Theorem 2.1.1 Fix an integer n0 > 0, then the Fourier coefficients of the volatility
are given by the following formulae:
N
π 1 2
a0 (Σ) = lim
N →∞ N + 1 − n0 ∑ ( a (dp j ) + bs2 (dp j )),
2 s
s = n0
N
2π
aq (Σ) = lim
N →∞ N + 1 − n0 ∑ (as (dp j )as+q (dp j )), (2.1.5)
s = n0
N
2π
bq (Σ) = lim
N → ∞ N + 1 − n0
∑ (as (dp j )bs+q (dp j )),
s = n0

∀q > 0 for aq (Σ) and ∀q ≥ 0 for bq (Σ).

By polarisation of the univariate case, the Fourier coefficients can be ex-
tended to the multivariate case, given by Theorem 3.2 in [1]:
Theorem 2.1.2 Fix an integer n0 > 0, then the Fourier coefficients of the volatility
are given by the following formulae:
N
π 1
a0 (Σi,j ) = lim
N →∞ N + 1 − n0 ∑ 2
( as (dpi ) as (dp j ) + bs (dpi )bs (dp j )),
s = n0
N
2π 1
aq (Σi,j ) = lim
N →∞ N + 1 − n0 ∑ 2
( as (dpi ) as+q (dp j ) + as (dp j ) as+q (dpi )),
s = n0
N
2π 1
bq (Σi,j ) = lim
N + 1 − n0 ∑ ( as (dpi )bs+q (dp j ) + as (dp j )bs+q (dpi )),
N →∞ s = n0 2
(2.1.6)
∀q > 0 for aq (Σij ) and ∀q ≥ 0 for bq (Σij ).
Remark 2.1.1 From deriving Theorem 2.1.1 and 2.1.2, we note that the scaling
factors are different to that of Theorem 3.1 and 3.2 in [1] but the same as [16].
Once the Fourier coefficients of the volatility matrix have been computed,
results from Fourier theory allows the reconstruction of Σ from its Fourier
coefficients [1]. Using the Fourier-Féjer inversion formula to reconstruct Σ,
we get
2.1. Trigonometric Fourier Transform 5

N  
k
Σ(t) = lim
N →∞
∑ 1−
N
( ak (Σ) cos(kt) + bk (Σ) sin(kt)) , ∀t ∈ (0, 2π ).
k =0
(2.1.7)
The Féjer inversion formula has the advantage that if Σ is a positive function,
then the approximation (2.1.7) will again be positive [1]. However, we are
more interested in the integrated volatility defined as

Z 2π
σ̂ij2 = Σi,j (t)dt. (2.1.8)
0

Which can be recovered by adjusting (2.1.5) to be

σ̂ij2 = 2πao (Σi,j ). (2.1.9)

2.1.2 Numerical Implementation

To implement this procedure, the first thing to do is to re-scale the irregularly
spaced observations [t1 , ..., tn ] into the interval [0, 2π ] [15] using the formula

2π (t j − t1 )
τj = , j = 1, ..., n. (2.1.10)
( t n − t1 )

The integrals for the Fourier coefficients of dp j can be computed using inte-
gration by parts [15]. Resulting in

p j (2π ) − p j (0) k 2π
Z
ak (dp j ) = + sin(kt) p j (t)dt,
π π 0 (2.1.11)
k 2π
Z
j
bk (dp ) = − cos(kt) p j (t)dt.
π 0

We note that (2.1.11) is numerically stable, because it does not involve the
differentiation of p j [1]. Furthermore, since the data gathered from financial
markets are discrete and therefore finite, we need an assumption of how the
data points are connected in order to compute (2.1.11). Malliavin and Man-
cino assume that p j (t) is equal to p j (ti ) in the interval [ti , ti+1 ] [1], also known
as the previous-tick interpolation [12]. Resulting in

N −1
p j (2π ) − p j (0) 1
ak (dp j ) ≈
π
+
π ∑ [cos(kti ) − cos(kti+1 )] p(ti ),
i =1
(2.1.12)
N −1
1
j
bk (dp ) ≈
π ∑ [sin(kti ) − sin(kti+1 )] p(ti ).
i =1
6 Chapter 2. The Theory

We further note that the choice of interpolation is important. Barucci and

Reno [7] showed that linear interpolation of prices between the interval [ti , ti+1 ]
resulted in a downward bias in the estimator. Malherbe further points out
this is because linear interpolation induces spurious auto-correlation [15].

Algorithm 1 Trigonometric Fourier Transform

Require:
1. (n x m) matrix P of asynchronously sampled price
2. (n x m) matrix T of asynchronously sampled times
Re-scale the time [tmin , tmax ] → [0, 2π ]
I. Extract trading times and prices
for i = 1 to m do
I.1. Slice the non-uniformly re-scaled sampled times
for the ith object
τ̃ ← τ (i )
I.2. Slice the sampled data indexing the times
for the ith object
ϕ̃ ← ln( p(τ̃ ))
II. Compute Fourier coefficients for all values of k
ϕ̃(2π )− ϕ̃(0) −1
ã ← π + π1 ∑ N
j=1 [cos( k̃ τ̃j ) − cos( k̃ τ̃j+1 )] ϕ̃ ( τ̃j )
−1
b̃ ← ∑N1
πj=1 [sin( k̃ τ̃j ) − sin( k̃ τ̃j+1 )] ϕ̃ ( τ̃j )
I.3. Gather Fourier coefficients for ith object
for all values values of k
a(i ) ← ã
b(i ) ← b̃
end for
Compute the integrated volatility and co-volatility over
the time window for objects i and j
2
Σij ← |πK | ∑k∈K [ ak (i ) ak ( j) + bk (i )bk ( j)]
Σij
Rij ← √ √
Σii Σ jj
return (Σ, R)

Algorithm 1 1 was provided by [17], [18] and [15].

2.2 Complex Exponential Fourier Transform

2.2.1 Derivation
We present the second derivation proposed by Malliavin and Mancino [2],
which uses a different Fourier approach to solve the same problem. The
assumptions required for the derivation is that p j (t) is a continuous semi-
martingale satisfying the stochastic differential equation
1 The pair-wise implementation can be found in ftcorr.R.
2.2. Complex Exponential Fourier Transform 7

d
∑ σi dW i + b j dt,
j j
dp = j = 1, ..., n, (A-I)
i =1

where W = (W 1 , ..., W d ) are independent Brownian motions on a filtered

probability space, σ∗∗ and b∗ are adapted stochastic processes satisfying

hZ T i
E (b j (t))2 dt < ∞,
0
hZ T i (A-II)
j
E (σi (t))4 dt < ∞,
0
i = 1, ..., d; j = 1, ..., n.

The idea main idea behind this derivation is the same as section 2.1. We
want to establish a connection between the Fourier transform of the volatility
process (2.1.2) and the Fourier transform of the price process [9].
We first re-scale the time window from [0, T ] to [0, 2π ]. We then define the
Fourier transform of dp j as

1
Z
j
F (dp )(k) := exp(−ikt)dp j (t), (2.2.1)
2π [0,2π ]

and the Bohr convolution product between two functions Φ, Ψ on the inte-
gers Z as

N
1
(Φ ∗ B Ψ)(k) := lim ∑ Φ ( s ) Ψ ( k − s ).
N →∞ 2N + 1 s=− N
(2.2.2)

This leads to Theorem 2.1 in [2]:

Theorem 2.2.1 Consider a process p satisfying the assumption (A-II). Then we
have for i,j = 1,2:

1
F (Σij )(k) = F (dpi ) ∗ B F (dp j )(k), ∀k ∈ Z. (2.2.3)
2π
The equality (2.2.3) is attained in probability, which means the limit in the convolu-
tion product exists in probability
Using Theorem 2.2.1, we get that

2π
F (Σij )(k) = lim
N →∞ 2N + 1
∑ F (dpi )(s)F (dp j )(k − s),
|s|≤ N (2.2.4)
∀k ∈ Z.
8 Chapter 2. The Theory

Now that we have an expression for the Fourier coefficients of the volatility
process, we can reconstruct Σ(t) using the Féjer inversion formula. Yielding

 |k| 
Σij (t) = lim
N →∞
∑ 1−
N+1
F (Σij )(k) exp(ikt). (2.2.5)
|k|≤ N

Having an expression for the Fourier coefficients of the volatility process al-
lows for the computation of the integrated volatility as

Z 2π
σ̂ij2 = Σij (t)dt = 2π F (Σij )(0)
0 (2.2.6)
2 i j
= (2π ) (F (dp ) ∗ B F (dp ))(0).

2.2.2 Numerical Implementation

We first re-scale the irregularly spaced observations [t1 , ..., tn ] to the interval
[0, 2π ] using (2.1.10). We then require an interpolation method for the dis-
crete observations for the price process. Malliavin and Mancino [2] use the
previous-tick interpolation and get

n −1
p1n (t) := ∑ p1 (t1i ) I[t1i ,t1i+1 ) (t),
i =1
n −1
(2.2.7)
p2n (t) := ∑p 2
(t2j ) I[t2 ,t2 ) (t).
j j +1
j =1

Malliavin and Mancino [2] defines Ii1 := [t1i , t1i+1 ) and Jj1 := [t2j , t2j+1 ) and the
returns by δI 1 ( p1 ) := p1 (t1i+1 ) − p1 (t1i ) and δ J 2 ( p2 ) := p2 (t2j+1 ) − p2 (t2j ). Then
i j
the Fourier coefficients of the price process through use of a simple function
approximation becomes

n −1
1
F (dp1n )(k) ≈
2π ∑ exp(−ikt1i )δIi1 ( p1 ),
i =1
n −1
(2.2.8)
1
F (dp2n )(k) ≈
2π ∑ exp(−ikt2j )δIj2 ( p2 ).
j =1

By combining (2.2.2), (2.2.6) and (2.2.8) together, we get

Z 2π n −1 n −1
1 is(t1 −t2 )
0
ij
Σ (t)dt = ∑ ∑ ∑
2N + 1 |s|≤ N j=1 i=1
e i j δ I 1 ( p1 ) δ I 2 ( p2 )
i j
(2.2.9)

for large n,N.

2.2. Complex Exponential Fourier Transform 9

Algorithm 2 Complex Exponential Fourier Transform

Require:
1. (n x m) matrix P of asynchronously sampled price
2. (n x m) matrix T of asynchronously sampled times
Re-scale the time [tmin , tmax ] → [0, 2π ]
I. Extract trading times and prices
for i = 1 to m do
I.1. Slice the non-uniformly re-scaled sampled times
for the ith object
τ̃ ← τ (i )
I.2. Slice the sampled data indexing the times
for the ith object
ϕ̃ ← ln( p(τ̃ ))
I.3. Compute the returns
δj ← ϕ̃(τ̃j+1 ) − ϕ̃(τ̃j )
II. Compute Fourier coefficients for all values of k
n−1 ikτ̃j
c̃+
k ← ∑ j =1 e δj
n−1 −ik τ̃j
c̃−
k ← ∑ j =1 e δj
I.3. Gather Fourier coefficients for ith object
for all values values of k
c+ (i ) ← c̃+
c− (i ) ← c̃−
end for
Compute the integrated volatility and co-volatility over
the time window for objects i and j
−
Σii ← |K1 | ∑k∈K [c+
k (i ) ck (i )]
−
Σij ← |K1 | ∑k∈K [c+
k (i ) ck ( j )]
Σ ji ← Σij
Σij
Rij ← √ √
Σii Σ jj
return (Σ, R)

Algorithm 2 2 was provided by [17] and [18].

Remark 2.2.1 The denominator of the scaling factors for a0 (Σ) and F (Σ)(0) de-
pends on how many Fourier coefficients we are summing over. For the Trig Fourier
transform, we sum from ∑sN=n0 therefore has N + 1 − n0 while the Fourier trans-
form, we sum from ∑|s|≤ N therefore has 2N + 1.
2 The pair-wise implementation can be found in ftcorr.R.
10 Chapter 2. The Theory

2.3 Hayashi-Yoshida Estimator

2.3.1 Derivation
We present the estimator proposed by Hayashi and Yoshida [3]. This estima-
tor has the ability to deal with asynchronous observations without the need
to synchronize the data beforehand. The estimator is an adjustment on the
well known realized covariance estimator defined as

n −1
V := ∑ ( Pt1i+1 − Pt1i )( Pt2i+1 − Pt2i ). (2.3.1)
i =1

The realized covariance estimator has the well known property that as the
RT
Mesh tends towards 0 (i.e. max1≤i≤n−1 |ti+1 − ti | → 0), then V → 0 Σij (t)dt
in probability. As Hayashi and Yoshida point out [3], there are two crucial is-
sues regarding the implementation of the realized covariance estimator. The
first is that actual transaction data is asynchronous. Secondly, due to the
asynchrony, a significant portion of the original data set will be missing at
pre-specified grid points. Therefore, in order to use (2.3.1), we must choose a
common interval h first, and impute or interpolate the missing observations
in some way [3] - this is commonly referred to as synchronizing the data. The
first thing to notice is that the estimate V heavily depends on the value of h
we pick. Additionally, as we have mentioned before Barucci and Renò [7]
found that linear interpolation induces a bias. Malherbe [15] points out that
the common intervals h need not be the same length, however we note that
it is important for the intervals to be common. Otherwise the contribution in
(2.3.1) will be 0.
Hayashi and Yoshida proposed a cumulative covariance estimator which is
free from the need to synchronize the data beforehand. The assumptions
required are that the price process follows the one-dimensional Itô process

dpl = µl dt + σl dW l , l = 1, 2, (A-III)

with dhW 1 , W 2 it = ρdt where ρ ∈ (−1, 1) is an unknown deterministic func-

tion, pl (0) > 0 is a constant, µl is a progressively measurable function and
σl > 0 is a deterministic and bounded function [3]. Furthermore, for the
sampling times. Let T ∈ (0, ∞) be an arbitrary terminal time for observing
the price processes. Let Π1 := ( I i )i=1,2,... and Π2 := ( J j ) j=1,2,... be the sets
of random intervals which partition (0, T ) for price process 1 and 2 respec-
tively. Let T 1,i := inf{t ∈ I i+1 } represent the ith observation time of P1 and
T 2,j := inf{t ∈ J j+1 } be that of P2 . Let n be the size of Π1 and Π2 . We assume
that the sampling intervals Π := (Π1 , Π2 ) satisfy the following

(i ). ( I i ) and ( J i ) are independent o f P1 and P2 .

(A-IV)
(ii ). As n → ∞, max | I i | ∨ max | J j | → 0.
i j
2.3. Hayashi-Yoshida Estimator 11

where | I | is the length of an interval I [3]. Hayashi and Yoshida define the
cumulative covariance estimator as

n n
Un := ∑ ∑ ∆P1 ( I i )∆P2 ( J j )1{ I i ∩ J j 6=∅} . (2.3.2)
i =1 j =1

This leads to Theorem 3.1 in [3]:

Theorem 2.3.1 Suppose assumption (A-IV) holds
(a) If sup0≤t≤T |µlt | ∈ L4 , l = 1, 2, then Un → θ as n → ∞.
(b) If sup0≤t≤T |µlt | < ∞ almost surely. l = 1, 2, then Un is consistent for θ, that is,
Un → θ in probability as n → ∞.
RT
Where θ := 0 σ1 σ2 ρdt = h P1 , P2 i T .
Remark 2.3.1 We note that (2.3.2) has n. However the size of Πi and Π j need not
be equal. We could have set Ni and Nj to be the size of Πi and Π j respectively [16].

2.3.2 Numerical Implementation

For the case when i = j the computation simply becomes the quadratic vari-
ation [16]. Therefore,

Z T Ni

0
Σii (t)dt = ∑ [∆P1 ( I i )]2. (2.3.3)
i =1

For the case when i 6= j, we use Kanatani’s weighted realized volatility [19]
defined as

Z T
0
Σij (t)dt = ∆Pi W∆P j
0
Ni Nj (2.3.4)
= ∑ ∑ wkl ∆P ( I i k j
)∆P ( J ),l

k =1 l =1

where

Pi (t1i ) − Pi (t0i )
   
w11 . . . w1Nj
..  . .. ..  .
∆Pi =   , W =  ..
 
. . . 
i i i i
P (t Ni ) − P (t Ni −1 ) w Ni 1 . . . w Ni Nj

The weights for Hayashi Yoshida are given by

(
j j
1 i f (tik−1 , tik ] ∩ (tl −1 , tl ] 6= ∅
wkl = (2.3.5)
0 otherwise.
12 Chapter 2. The Theory

Remark 2.3.2 Kanatani’s weighted realized volatility can also be used for Malliavin
and Mancino’s Fourier estimator. See [19] or [16].

Algorithm 3 Hayashi Yoshida

Require:
1. (n x m) matrix P of asynchronously sampled price
2. (n x m) matrix T of asynchronously sampled times
Loop through every element of Σ [m x m]
for i = 1 to m do
I.1. Compute the returns
δi ← ln( pi (tk )) − ln( pi (tk−1 ))
for j = 1 to m do
δj ← ln( p j (tk )) − ln( p j (tk−1 ))
I.2. Compute Kanatani’s weight matrix for the
i-th and j-th stock
W ← Kanatani weight for HY
I.3. Compute Σij
Σij ← δi0 Wδj
end for
end for
II. Compute the correlation matrix
Σij
Rij ← √ √
Σii Σ jj
return (Σ, R)

Algorithm 3 3 was provided by [18].

3 The pair-wise implementation can be found in ftcorr.R.

 13

Chapter 3

Monte Carlo Experiments

Monte Carlo experiments are conducted to identify the differences between

the two estimator in hope to determine which estimator is better between the
two.

3.1 Effect of missing observations

The first experiment focuses on how the two estimators differ when asyn-
chronicity is induced by down-sampling the price path. The experiment is
conducted by simulating 10,000 seconds from a bivariate Geometric Brown-
ian motion with daily parameters µ1 = 0.01, µ2 = 0.01, σ12 = 0.1, σ22 = 0.2, ρ12
ranging from (−1, 1) and a starting price of R100. The non-synchronicity is
achieved by randomly sampling a percentage of the observations from each
sample path and removing them. The Geometric Brownian motion satisfies
the following system of SDEs

dSi (t)
= µi dt + σi dWi (t), i = 1, 2. (3.1.1)
Si ( t )

Figure 3.1 1 (a), we see that both MM (blue dotted line) and HY (red dotted
line) perfectly recover the induced correlation (black dotted line) for the syn-
chronous case. From figure 3.1 (b) through to (d), it is clear that as the level of
asynchrony increases, MM appears to have a downward bias towards zero
which [12] attributes to the Epps effect [4] while HY recovers the induced
correlation regardless of the level of asynchrony.
Hayashi and Yoshida claim that the Epps effect is a bias that arises from the
estimator for which their estimator is immune to [3]. Looking at figure 3.1,
this seems to be the case. However, this goes against the findings of [6],
[12], [13], [20]. The current literature has identified the main sources for the
Epps Effect to be: smaller sampling intervals [4], [20], lead-lag [6], [12] and
asynchronicity [12], [13]. Closed-form expressions recovering the Epps Effect
can be found in [6], [20] - indicating that the Epps Effect is not a bias from the
1 Figure 3.1 can be reproduced using MissingData.R.
14 Chapter 3. Monte Carlo Experiments

estimator. This, in turn, means that it is the HY that is upward biased even
though it recovers the induced correlation.

(a) Synchronous (b) Asynchronous (10%) Missing

1.0 1.0

0.5 Method 0.5 Method

correlation (ρ)

correlation (ρ)
Induced Induced

0.0 0.0
MM MM

HY HY
−0.5 −0.5

−1.0 −1.0
0 5 10 15 20 0 5 10 15 20
simulation simulation

(c) Asynchronous (20%) Missing (d) Asynchronous (40%) Missing

1.0 1.0

0.5 Method 0.5 Method

correlation (ρ)

Induced correlation (ρ) Induced

0.0 0.0
MM MM

HY HY
−0.5 −0.5

−1.0 −1.0
0 5 10 15 20 0 5 10 15 20
simulation simulation

F IGURE 3.1: Comparing different levels of missing data to

demonstrate the bias between the MM and HY estimators. Con-
cretely, (a) through to (d) show 0%, 10%, 20% and 40% of each
sample path replaced by missing data to replicate asynchrony.
As per the figure legend, the blue dotted line is the MM esti-
mator, estimated using Algorithm 2; the red dotted line is the
HY estimator, estimated using Algorithm 3 and the black dot-
ted line is the induced correlation of the GBM. Each simulation
with varying correlation is done by simulating 10,000 seconds
from a bivariate GBM satisfying (3.1.1) and numerically sim-
ulated using Algorithm 10. The figure shows that if the data
is discrete and asynchronous then the MM estimator is most
appropriate, however, if the data is sampled discretely from
an underlying continuous-time GBM then the HY estimator is
more appropriate. We argue that the Epps Effect from asyn-
chrony is real and as such due to the discrete nature of financial
market trades one should use the MM estimator rather than the
HY estimator.

This experiment although recovers the Epps effect arising from asynchrony,
is not an experiment conducted on a truly asynchronous process. It is rather
a missing observation experiment. Therefore, the argument is that the HY
estimator is the better estimator of the two if one believes that the observed
prices in the market are discrete samples of an underlying continuous stochas-
tic process and that asynchrony is a missing data problem. Then the HY
estimator will be able to reproduce the true underlying correlation between
3.2. Effect of the SDE 15

the assets by allowing multiple contributions to the estimator, which would

imply the MM estimator has a downward bias attached to it. On the other
hand, the argument is that the MM estimator is the better estimator of the
two if one is of the belief that the world is not a missing data problem from an
underlying synchronous continuous stochastic process; but is rather discon-
nected, discrete and asynchronous where the events and their relationships
to each other are the fundamental entities and measurables of the finance
world. Then the MM estimator will produce the true correlation in the sys-
tem as it is lossless interpolation between the events. This would imply that
the HY estimator will have an upward bias that is caused by the multiple
contributions [17].
We argue the Epps effect is a fundamental property of financial market data
which is not picked up by the HY estimator, and furthermore, TAQ data is
truly discontinuous [21], discrete and asynchronous events which although
does not fit into the framework of the two estimators; MM is the best tool we
have for studying co-movements between discrete events due to the lossless
interpolation. Therefore we argue MM is the preferred estimator of the two
when it comes to studying high-frequency data.
Additional issues regarding the HY estimator is pointed out in [10], [11]. The
first issue with the HY estimator is that when the processes are highly asyn-
chronous, the HY estimator deletes observations through its multiple contri-
butions (e.g. Fig. 1 in [11]), therefore it does not utilise all available obser-
vations. Furthermore, a critical assumption underlying the HY estimator is
that the correlation between two assets does not extend beyond the intervals
where returns fully or partially overlap. Meaning that information regarding
the correlation is fully accounted for when a price update arrives. This as-
sumption does not hold in practice and causes the HY estimator to be biased
as shown by [10].

3.2 Effect of the SDE

The first experiment found that the two estimators differ under asynchrony,
and the level of asynchrony determines how different the estimators behave.
To further gain insight into the two estimators, various stochastic processes
are studied to identify alternative situations where these two estimators dif-
fer.
The second experiment focuses on if alternative stochastic processes will
cause the two estimators to differ. To this end, the Merton model, Variance
Gamma, GARCH (1,1) and Ornstein Uhlenbeck will be used to compare the
two estimators.
The bivariate Merton model satisfies the following system of SDEs

dSi (t)
= µi dt + σi dWi (t) + dJi (t), i = 1, 2. (3.2.1)
Si (t−)
16 Chapter 3. Monte Carlo Experiments

where corr(dW1 , dW2 ) = ρ. The Ji are independent of the Wi with piece-wise

constant paths [22]. J is defined as

(a) λ = 0 (b) λ = 0.2

1.0 1.0
Method Method

0.5 Induced 0.5 Induced

correlation (ρ)

correlation (ρ)
MMSyn MMSyn

0.0 0.0
HYSyn HYSyn

MMAsyn MMAsyn
−0.5 −0.5
HYAsyn HYAsyn

−1.0 −1.0
0 5 10 15 20 0 5 10 15 20
simulation simulation

(c) λ = 0.5 (d) Variance Gamma

1.0 1.0
Method Method

0.5 Induced 0.5 Induced

correlation (ρ)
correlation (ρ)

MMSyn MMSyn

0.0 0.0
HYSyn HYSyn

MMAsyn MMAsyn
−0.5 −0.5
HYAsyn HYAsyn

−1.0 −1.0
0 5 10 15 20 0 5 10 15 20
simulation simulation

F IGURE 3.2: Comparing various stochastic processes to demon-

strate the effect on the MM and HY estimators. Concretely,
(a) shows a pure diffusion process, (b) and (c) show a jump-
diffusion process and (d) shows a pure jump process. As per
the figure legend, the blue and purple dotted lines are the MM
estimator estimated using Algorithm 2 under synchronous and
asynchronous observations respectively. The red and orange
dotted line is the HY estimator estimated using Algorithm 3 un-
der synchronous and asynchronous observations respectively.
The black dotted line is the induced correlation from the var-
ious stochastic processes. (a) through to (c) is 10,000 seconds
simulated from a bivariate Merton model satisfying (3.2.1) and
numerically simulated using Algorithm 11. (d) is from a bi-
variate Variance Gamma model satisfying (3.2.3) and numeri-
cally simulated using Algorithm 13. The asynchrony is induced
by down-sampling 20% of the observations from each sample
path. The figure shows that both MM and HY produce the same
estimates regardless of the underlying process and that it seems
the difference between them arises due to asynchrony. Under
asynchrony, HY recovers the synchronous estimates while the
correlation for MM drops.

N (t)
Ji (t) = ∑ (Yj − 1), (3.2.2)
j =1
3.2. Effect of the SDE 17

where N (t) is a Poisson process with Yj ∼ LN ( a, b) i.i.d and also independent

of N (t).
The daily parameters used for the Merton model are µ1 = 0.01, µ2 = 0.01,
σ12 = 0.1, σ22 = 0.2, ρ12 ranging from (−1, 1), a1 = 0, a2 = 0, b1 = 100,
b2 = 100 and varying λ to move from a pure diffusion process to a jump-
diffusion process. A sample path of 10,000 seconds is simulated starting at
R100. This model will determine if the two estimators differ due to the effect
of jumps.
The bivariate Variance Gamma (VG) process satisfies the following SDEs:

Si ( t ) = U ( t ) − D ( t ) , i = 1, 2, (3.2.3)

with U and D being independent gamma processes satisfying

U (ti+1 ) − U (ti ) ∼ Gamma(α(ti+1 − ti ), β),

(3.2.4)
D (ti+1 ) − D (ti ) ∼ Gamma(α(ti+1 − ti ), β).

U and D are limited to have the same shape and scale parameters allowing
an alternative representation W ( G (t)) where W is a standard Brownian mo-
tion, G a gamma process [22] and corr(dW1 , dW2 ) = ρ ranging from (-1, 1).
A sample path of 10,000 seconds is simulated starting at R100 with daily pa-
rameters µ1 = µ2 = 0.01, σ12 = 0.1, σ22 = 0.2 and β 1 = β 2 = 1. This model
will determine if a pure jump process will cause the estimators to differ.
The bivariate GARCH (1,1) model satisfies the following SDEs:

dpi (t) = σi (t)dWi (t), i = 1, 2, (3.2.5)

and
p
dσ12 (t) = θ1 [w1 − σ12 ]dt + 2λ1 θ1 σ12 (t)dW3 (t),
p (3.2.6)
dσ22 (t) = θ2 [w2 − σ22 ]dt + 2λ2 θ2 σ22 (t)dW4 (t).

where corr(dW1 , dW2 ) = ρ ranges from (-1, 1). We simulate a sample path
of 10,000 seconds starting at R100 using the parameters from [12], [23] i.e.
θ1 = 0.035, θ2 = 0.054, w1 = 0.636, w2 = 0.476, λ1 = 0.296 and λ2 = 0.48
2 . This model will determine if stochastic volatility and volatility clustering

will cause the estimators to differ.

The bivariate Ornstein Uhlenbeck process satisfies the following SDEs:

dpi (t) = θi (µi − pi (t))dt + σi dWi (t), i = 1, 2, (3.2.7)

2 We note the SDE specified by [12] is different to what [23] has, we followed [23].
18 Chapter 3. Monte Carlo Experiments

where corr(dW1 , dW2 ) = ρ ranges from (-1, 1). A sample path of 10,000 sec-
onds is simulated starting at R100 with parameters µ1 = µ2 = 100, σ12 = 0.1,
σ22 = 0.2, θ1 = 0.035 and θ2 = 0.054. This model will see if the two estimators
differ under mean-reversion.
(a) GARCH (1,1) (b) Ornstein Uhlenbeck
1.0 1.0
Method Method

0.5 Induced 0.5 Induced

correlation (ρ)

correlation (ρ)
MMSyn MMSyn

0.0 0.0
HYSyn HYSyn

MMAsyn MMAsyn
−0.5 −0.5
HYAsyn HYAsyn

−1.0 −1.0
0 5 10 15 20 0 5 10 15 20
simulation simulation

F IGURE 3.3: Comparing the two estimators based on volatil-

ity clustering (a) and mean-reversion (b). As per the figure
legend, the blue and purple dotted lines are the MM estima-
tor estimated using Algorithm 2 under synchronous and asyn-
chronous observations respectively. The red and orange dot-
ted line is the HY estimator estimated using Algorithm 3 un-
der synchronous and asynchronous observations respectively.
The black dotted line is the induced correlation from the vari-
ous stochastic processes. (a) is 10,000 seconds simulated from a
bivariate Geometric GARCH(1,1) satisfying (3.2.5) and numer-
ically simulated using Algorithm 12. (b) is from a bivariate Ge-
ometric Ornstein Uhlenbeck satisfying (3.2.7) and numerically
simulated using Algorithm 14. The figure shows that both es-
timators produce the same estimates in the synchronous case
while for the asynchronous case, HY recovers the synchronous
estimates but the correlation drops for MM.

Figure 3.2 3 (a) through to (c), λ1 = λ2 = 0, 0.2 and 0.5 respectively. The
asynchronicity is induced by down-sampling each price path by 20%. For all
the plots in figure 3.2, the synchronous MM (blue dotted line) is the same as
the synchronous HY (red dotted line). The asynchronous HY (orange dotted
line) recovers the synchronous estimates while the asynchronous MM (pur-
ple dotted line) has a lower correlation estimate than the synchronous case.
In figure 3.2 (a) through to (c), both the synchronous MM and HY estimators
produce the same estimate which drops towards zero as λ increases. This
drop in correlation is not due to any bias from the estimators but rather due
to the fact that the jump process of the Merton model is independent of the
underlying diffusion process. Therefore as the intensity of jumps increase,
the impact from the independence, seeps through to change the correlation
structure of the overall jump-diffusion process. This is the case because when
the trades are synchronous, HY becomes the Realized Volatility (RV) [16] and
the RV is consistent under jumps [24]. Therefore, the two estimators seem to
3 Figure 3.2 can be reproduced using SDE1.R.
3.3. Effect of Asynchrony 19

only differ under asynchrony and it is not dependent on the type of diffusion
process.
Figure 3.3 4 , the asynchrony is induced by down-sampling each price path
by 20%. Once again, the synchronous MM (blue dotted line) is the same
as the synchronous HY (red dotted line) and the asynchronous HY (orange
dotted line) recovers the synchronous estimates while the asynchronous MM
(purple dotted line) has a lower correlation estimate than the synchronous
case. Figure 3.3 provides the insight that neither volatility clustering nor
mean-reversion causes the two estimators to differ under synchronicity. The
two estimators only seem to produce different estimates under asynchronous
conditions.
This experiment falsifies the idea that various stochastic processes will cause
the two estimators to differ, rather it further validates that the two estima-
tors differ only under asynchronous conditions, where the HY estimator is
immune to the Epps effect brought through by asynchrony (under missing
data conditions) while MM estimator picks up the Epps effect.

3.3 Effect of Asynchrony

In the previous experiments, asynchrony is induced by removing observa-
tions from the sample path. Although asynchrony is achieved; it is more of a
missing data problem. The next experiment will follow a similar methodol-
ogy used by [12] and [13]. The focus will be to achieve asynchronous sample
paths that behave more like tick-by-tick TAQ data and investigate the effect
the number of Fourier coefficients (N in (2.2.2)) has on the estimates.
Naturally, we first recover the results from [12]. The first thing to point out
is that [12] has a different specification of the GARCH (1,1) compared to [23]
- where the parameters were borrowed from. In figure 3.4 we present both
specifications of the GARCH (1,1). The experiment in figure 3.4 is conducted
by first simulating price paths of 86,400 seconds from a bivariate GARCH
(1,1) with parameters from above. The asynchrony is induced by sampling
the price path with an exponential inter-arrival time with a mean of 15 sec-
onds and 45 seconds from asset 1 and asset 2 respectively. The synchronous
case here is achieved by forcing the first time series to be observed at the
same time as the second time series (i.e. the price paths were sampled with
the same exponential inter-arrival time with a mean of 45 seconds). For each
of the asynchronous and synchronous cases, we compute the correlation esti-
mate using algorithm 2 for N ranging from 10 to 160. The effect achieved by
studying the ranging Fourier coefficients (N) is that it allows various sam-
pling frequencies to be studied through the relationship of the two from
Fourier analysis.
4 Figure 3.3 can be reproduced using SDE2.R.
20 Chapter 3. Monte Carlo Experiments

From figure 3.4 5 , the Epps effect is clearly demonstrated. As N increases, so

does the sampling frequency and from the asynchronous case (green dots) it
is clear that the correlation drops as the sampling frequency increases. How-
ever, for the synchronous case (blue dots), the correlation does not drop as
the sampling frequency increases. This is because the stochastic processes we
have studied have dimensionless correlation (independent of the sampling
intervals). The Epps effect recovered here is due to asynchrony which is dif-
ferent from the original Epps effect presented by Thomas Epps [4], where the
correlation drops from synchronous observations as the sampling interval
decreases. Researchers have investigated this - specifically [20] was able to
derive an analytical expression for the Epps effect arising from smaller sam-
pling intervals by decomposing the correlation of time scale ∆t as a function
of lagged autocorrelations and correlations of smaller time scales ∆t0 . Ad-
ditionally, [6] was able to further extend the results from [12] by analytically
deriving the Epps effect arising from asynchrony as a function of ∆t.

(a) GARCH (1,1) − Reno (b) GARCH (1,1) − Andersen

0.35
0.355

Method Method
0.350
Correlation (ρ)

Correlation (ρ)

Induced Induced
0.34

0.345
MM Asyn MM Asyn

MM Syn
0.340 MM Syn
0.33

0.335

40 80 120 160 40 80 120 160

Max Fourier Coefficients (N) Max Fourier Coefficients (N)

F IGURE 3.4: We recover the result from [12] using the complex
exponential Fourier Transform 2. The average correlation as
a function of the sampling frequency N in (2.2.2). Concretely,
the asynchronous sample paths for (a) and (b) are exponen-
tial inter-arrival time samples from 86,400 seconds of simulated
data. The exponential inter-arrival times have a mean of 15 sec-
onds and 45 seconds respectively for asset 1 and asset 2. The
synchronous sample paths for (a) and (b) are achieved by forc-
ing the first time series to be observed at the same times as the
second time series. (a) is simulated by adjusting (3.2.5) to how
[12] defined the SDE and implemented by adjusting algorithm
12 accordingly. (b) is simulated from (3.2.5) using algorithm 12.
As per the figure legend, the green dots and blue dots are the
asynchronous and synchronous sample paths estimated using
algorithm 2 respectively. The orange line is the induced corre-
lation between (3.2.5). The results are obtained through 10,000
replications.

Drawing attention back to comparing the two estimators, we modify the ex-
periment slightly. Specifically, for figure 3.5, the experiment is conducted
5 Figure 3.4 can be reproduced using Reno Recovery.R.
3.3. Effect of Asynchrony 21

by first simulating price paths of 10,000 seconds from the various stochas-
tic processes from above, using their respective parameters as before 6 . The
asynchrony is induced by sampling the first asset with an exponential inter-
arrival time with mean 30 seconds and the second asset with a mean of 45
seconds. The synchronous case here is achieved by forcing the first time se-
ries to be observed at the same time as the second time series. The rationale
behind adjusting the experiment is so that the Nyquist frequency can be in-
dicated.
From figure 3.5 7 , it is clear that for the MM estimates, the correlations de-
crease for the asynchronous case as the number of Fourier coefficients (N)
increase; whereas for the synchronous case the correlations become closer
to the synchronous HY estimates as N increases. Additionally, the error
bars calculated to be the standard deviation from the estimates decrease as
N increases indicating that the estimates become more accurate with more
Fourier coefficients. The HY estimates are not a function of N, but rather it
a baseline to compare the MM estimate against. For figure 3.5 (a) through
to (e), the asynchronous HY estimate recovers the synchronous HY estimate
which is expected as Hayashi and Yoshida have claimed that their estimator
is immune to the Epps effect. Oddly enough, the HY estimator demonstrates
an Epps effect when using the Ornstein Uhlenbeck process. This was not
picked up by the experiments before and a possible explanation for this is
due to how this experiment is set up. In the previous experiments asyn-
chrony was induced through a missing data manner, while in this experi-
ment, the asynchrony is induced through exponential inter-arrival times to
sample the price paths. Combined with the mean-reversion from the OU
process, the sampling may have picked up different co-movements between
the price paths. Another possible explanation is due to the combination of
sampling method and mean-reversion, spurious lead-lag relations may have
arisen due to the high levels of asynchrony which is another source for the
Epps effect as investigated by [6], [12]. This further highlights the downfall of
the HY estimator for high-frequency finance. Although the HY estimate may
be immune to the Epps effect in a missing data manner, it is not immune to
the Epps effect arising from lead-lag [10] nor from smaller synchronous sam-
pling intervals [4] as seen in real financial data in later sections. Additionally,
when levels of asynchrony is high - common for high-frequency data, the
HY deletes observations [11] and therefore it is not well suited to study the
co-movement between events.
6 The Merton Model has λ1 = λ2 = 0.2
7 Figure 3.5 can be reproduced using Reno Extended.R.
22 Chapter 3. Monte Carlo Experiments

(a) GARCH (1,1) − Reno (b) GARCH (1,1) − Anderson

Method Method
0.5 0.5
HY Asyn HY Asyn

0.4 0.4
Correlation (ρ)

Correlation (ρ)
HY Syn HY Syn

0.3 0.3
Induced Induced

0.2 0.2
MM Asyn MM Asyn
Nyquist Frequency Nyquist Frequency

0.1 0.1
MM Syn MM Syn

40 80 120 160 40 80 120 160

Max Fourier Coefficients (N) Max Fourier Coefficients (N)

(c) GBM (d) Merton Model

Method 0.5
Method
0.5
HY Asyn HY Asyn

0.4
0.4
Correlation (ρ)

Correlation (ρ)
HY Syn HY Syn

0.3
0.3
Induced Induced

0.2
0.2
MM Asyn MM Asyn
Nyquist Frequency
0.1 Nyquist Frequency
0.1
MM Syn MM Syn

40 80 120 160 40 80 120 160

Max Fourier Coefficients (N) Max Fourier Coefficients (N)

(e) Variance Gamma (f) Ornstein Uhlenbeck

Method Method
0.6

HY Asyn HY Asyn

0.50
Correlation (ρ)

Correlation (ρ)

HY Syn HY Syn
0.4

Induced Induced
0.25

MM Asyn MM Asyn
0.2 Nyquist Frequency

0.00 Nyquist Frequency

MM Syn MM Syn

40 80 120 160 40 80 120 160

Max Fourier Coefficients (N) Max Fourier Coefficients (N)

F IGURE 3.5: Comparing the two estimators by introducing

asynchrony through sampling the time series with exponen-
tial inter-arrival times. The parameters used are the same as
the previous experiments. We note that (a) uses the specifica-
tion of GARCH (1,1) by [12] while (b) uses the specification of
[23]. As per the legend, the synchronous (pink dots) and asyn-
chronous (blue dots) average MM estimates using algorithm 2
differ more as the number of Fourier coefficients (N) increase.
The synchronous (dark green line) and the asynchronous (or-
ange line) average HY estimates using algorithm 3 in general
recover the same estimates except for (f) - the OU process. The
error bars are plotted as the standard deviation of the estimates
across the replications. The green line is the induced correlation
for each of the processes set to be ρ = 0.35. Additionally, the
Nyquist frequency is indicated by the black dashed line. The
results are obtained through 1,000 replications.

From figure 3.5, the Nyquist frequency is calculated using the average sam-
pling frequency, this is not the true cutoff required to avoid any aliasing.
3.3. Effect of Asynchrony 23

The true cutoff required is computed by first finding the highest sampling
frequency present in the data, then computing the corresponding Nyquist
frequency. Using the true cutoff is how we compute all the MM estimates in
this paper except in figure 3.4 and 3.5. The rationale behind this is that we
are trying to study the co-movement between high-frequency events, there-
fore picking a lower N to avoid market microstructure noise will lead to the
aliasing of the event data. Picking a lower N, in essence, creates a smooth-
ing effect due to aliasing of higher frequencies which is useful to identify the
true signal under the market microstructure noise argument [9], but for the
context of identifying the co-movement of events, picking a lower N will be
a fatal choice to make.
This means that using the appropriate cutoff for the asynchronous cases in
figure 3.5, all the correlations for the MM estimate diminish to zero. This
result combined with the drop in correlation as the % of missing data in-
creased in figure 3.1, indicates the importance the level of asynchrony has in
contributing towards the Epps effect [13].
A point to be noted is that although this paper argues for the efficacy of
the MM estimator over the HY estimator in the high-frequency paradigm
through an event-based view of the world, we have presented results from
the classical continuous-time stochastic processes which is a slight inconsis-
tency to the event-based view we are taking. This is because of the limited
techniques present in the literature to simulate a price process. Therefore
we do note an extension on our results is to study how the correlations be-
have when the process is simulated from a Hawkes process [25] which can
hopefully provide more insight into the co-movement between events and
their relation to the Epps effect. However, our results are not futile because
we have recovered and validated the work from previous researchers and
have further performed a wholistic comparison of the two estimators with
the various SDEs.
Although we have argued for the efficacy of the MM estimator over the HY
estimator through an event-based view of the world, other researchers have
also argued for the efficacy of the MM estimator in the market-microstructure
noise view of the world. Specifically, [9] shows that the MM estimator is
unbiased for the contaminated price process by an appropriate choice of n
and N; while [11] points out that the HY estimator is infeasible in the setting
of market microstructure noise.
Finally, a subtle point to notice is that all the SDEs used in this paper has
dimensionless correlation which does not depend on time and therefore we
could not study the Epps effect arising from smaller sampling intervals us-
ing Monte Carlo experiments, but from figure 3.4 and 3.5, we use the Fourier
methods to study smaller sampling intervals. This subtle difference is due
to what [6] showed. The correlation arising from asynchrony not only de-
pends on the level of asynchrony, but also on the sampling intervals chosen.
Therefore it seems that the Epps effect arising from smaller sampling inter-
vals and asynchrony have some form of relation, and therefore more rigorous
24 Chapter 3. Monte Carlo Experiments

research into this aspect is required - by attempting to decompose these two

factors contributing towards the Epps effect. However, a possible complica-
tion that will arise is that [6] analytically showed the Epps effect arising from
asynchrony as a function of ∆t given by:
 
1  −λ∆t
ρ̃12
∆t = c 1+ e −1 , (3.3.1)
λ∆t

where the correlation only decreases and does not change signs. However,
in section 5, we show that there seems to be a structural change in the corre-
lation that it is not only decreasing but also becoming positively correlated.
Indicating that there is more to this problem than what meets the eye.
 25

Chapter 4

Data Engineering

Managing the data is a difficult problem which is often downplayed, there-

fore we begin by describing the challenges faced with managing high-frequency
data and how they were resolved. Furthermore, we outline all the algorithms
employed to create the datasets used in the analysis. All the aggregation
methods were built from the bottom up to ensure a white box process is
adapted for reproducible research.

4.1 Data Collection

4.1.1 Collecting Data From Bloomberg
The data obtained for the analysis consists of Trade and Quote (TAQ) data
from 10 equities listed on the Johannesburg Stock Exchange (JSE). The period
considered is the week from 31/05/2019 to 07/06/2019. The equities consid-
ered are FirstRand Limited (FSR), Shoprite Holdings Ltd (SHP), Absa Group
Ltd (ABG), Nedbank Group Ltd (NED), Standard Bank Group Ltd (SBK),
Sasol Ltd (SOL), Mondi Plc (MNP), Anglo American Plc (AGL), Naspers Ltd
(NPN) and British Am. Tobacco Plc (BTI). These equities are chosen due to
their high liquidity, because we are interested in high-frequency event data.
There exists various APIs one can use to collect data from Bloomberg. These
are outlined in table 4.1 provided by [26]. Each of the methods presented
in table 4.1 have their strengths and weaknesses depending on the type of
permission access one has on the shared Bloomberg terminal.

API Benefits Issues

Manual (GUI) Good for content discovery Documentation of user actions difficult.
Scripting though VBA is slow.
Good for human-directed data refresh.
Excel Add-In Only flat, tabular data or
Next-best for exploration
hard-to-machine-read data structures
Fast, can retrieve multidimensional data objects, Requires admin rights; low-level language.
C API
scriptable Poor for exploration.
Depends on C under hood;
Python API Easy language, can retrieve data objects, scriptable requires admin rights.
Poor for exploration.
R does not require admin rights,
Fast(er than Excel), easy language, mature IDE, can
R API but installing packages do.
retrieve data objects, scriptable
Poor for exploration.

TABLE 4.1: Bloomberg APIs compared

26 Chapter 4. Data Engineering

The manual (GUI) and excel add-in are unreliable when it comes to extract-
ing large datasets from Bloomberg; furthermore, it has the added complica-
tion of not being easily reproducible. Thus we will consider more dynamic
programmatic methods such as C, Python and R. The main complication that
arises with C and Python is that one might not have the administrative rights
to run the APIs, therefore we opt for the R API to extract TAQ data from
Bloomberg.
For the purpose of collecting TAQ data, the choice of R comes with its caveats,
mainly:
• Writing the TAQ data into flat files results in large flat files, which can
take several hours to complete depending on the memory available on
the terminal.
• Large flat files are complicated to read in R and often require Java or
C++ interfaces to speed up the process.
This issue is not trivial due to the large nature of TAQ data. For example,
one of the more liquid tickers - Naspers (NPN) has 15,544,244 data points
for a period of 6 months which results in a flat-file of 450-500MB. Therefore
obtaining data for all 10 tickers translates to roughly 5GB of data, illustrating
the non-trivial nature of this problem.
After identifying these issues, we realised the need to find a more efficient
and easily reproducible process to overcome the issues presented above.

4.1.2 Updated HF-Data Pipeline

The focus in our approach was to have an easily reproducible process that
spends minimal time on the Bloomberg terminal and can easily be extended.
In figure 4.1, we illustrate the data pipeline.

Bloomberg Load Data

R Environments
of Ticker Data

Data Processing

Data Pulling

F IGURE 4.1: Data Collection Pipeline

4.2. Data Cleaning 27

The key factor that solves the issue of large files and long computation time
is simply saving the data as R environments rather than flat files. The advan-
tages achieved by doing so are exceptional, namely:
• R environments are significantly smaller when saving data compared
to flat files. For example, the Naspers (NPN) ticker when saved to a
flat-file results in a size of 500MB, while as an R environment, the size
reduces to 47MB - a decrease of 90.6% in storage size.
• Reading in the R environment is significantly quicker - even more so
than using third party interfaces such as Java or C++ and furthermore
avoids any memory issues which arise from large file sizes.
Now the largest overhead left is simply the time it takes to extract data from
Bloomberg since the other areas have been streamlined. To see the efficiency
up the updated pipeline, we consider the file sizes and time spent in extract-
ing and loading TAQ data. We initially pulled 24 tickers which took approx-
imately 3 hours to extract from Bloomberg. Now saving the extracted data
as an R environment resulted in a file of 557.6MB as opposed to the 4.9GB
when saving the data as flat files. The point of significance is when we have
to read the data. Loading the R environment takes on average 25-30 seconds
while loading the flat files took several hours with no result 1 . Finally, the last
advantage this approach presents is that it is ring-fenced within R and does
not rely on third-party interfaces such as rJava or Rcpp 2 .
The pipeline although advantageous, does present some potential pitfalls
which we have not yet encountered. Specifically, if the R environment ex-
ceeds 2GB then reading in data might pose an issue. This can be pragmati-
cally solved by writing various assets into their own R environment. How-
ever, this should never present itself as a real issue given the limitation of
the Bloomberg terminal which only allows for 6 months of TAQ data to be
extracted. Finally, the pipeline is designed to extract TAQ data, but this can
easily be extended to pull other forms of data from Bloomberg.

4.2 Data Cleaning

4.2.1 Data Types
From the TAQ data, there are three types of observables: the bid, ask and
actual trades. We are only interested in the actual trades as they form the
price paths of interest. Furthermore, there are different trade types, exam-
ples include Automated Trade (AT), Late Trade (LT), Post Contra Trade (LC)
and Indicative Auction Information (IP) [27]. For the analysis, we extract
only the Automated Trades because only these form the continuous trading
1 The reason it did not finish loading is because the file size was too large, coupled with
the fact that R reads the entire dataset into RAM all at once, lack of RAM resulted in the
loading to never finish. Solving this requires the package bigmemory to overcome the issue
- the data gets read into the HDD, which comes with its own caveats.
2 bigmemory is part of the Rcpp family.
28 Chapter 4. Data Engineering

process. The other trade types are after-hour trades (LT), correction of previ-
ous days published off book trade (LC) and an indicative auction price based
on the volume maximising auction algorithm used to determine the auction
uncrossing price (IP) [27] - which are irrelevant to the analysis.

4.2.2 Aggregation
An issue with Bloomberg data which is not present with Thomson Reuters
data is that timestamps are only shown up to seconds, therefore there are
multiple trades with the same timestamp - illustrated in figure 4.2. This poses
an issue when using the MM and HY estimators - the two estimators require
unique time stamps for each trade.

F IGURE 4.2: Trades with the same time stamp.

To overcome this issue we first need to aggregate these “repeated” 3 trades.

The aggregation algorithm is presented in algorithm 4 4 and it uses a Volume
Weighted Average Price (VWAP) method of aggregation. VWAP was em-
ployed because it gives a better representation of the data given by the fact
that it weights each trade by the volume, which is directly linked to the price
impact [28].
3 The term repeated is used very loosely as these are not actually repeated trades.
4 The implementation can be found in AsynchronousData.R.
4.2. Data Cleaning 29

Algorithm 4 Aggregation of Repeated Trades

Require:
1. Ti the trading times, i = 1, ..., N
2. Si the observed prices, i = 1, ..., N
3. Vi the volume associated with the trade, i = 1, ..., N
Identify the unique trading times t∗j , j = 1, ...M
Gather trades with the same trading times into a set Jj , j = 1, ..., M
Procedure for the jth set:
1. Set
∑i∈ Jj pricei ∗ volumei
s∗j =
∑i volumei
2. Set Vj∗ = ∑i∈ Jj Vi
return ( T ∗ = {t∗j } jM=1 , S∗ = {s∗j } jM=1 , V ∗ = {Vj∗ } jM=1 )

In figure 4.3 we demonstrate the output of algorithm 4 using one set of “re-
peated” trades Jj from figure 4.2.

F IGURE 4.3: Trades aggregated into a unique time stamp.

4.2.3 Overnight returns

The trading times for the JSE begin at 09:00 and end at 17:00 5 , while the
closing auction begins at 16:50. Furthermore, the opening prices and clos-
ing prices are determined using a volume maximising auction algorithm
to determine the uncrossing price, therefore these prices not determined by
Market Orders (MO) hitting the Limit Order (LO) - which form the continu-
ous trading process, but rather these prices follow a Walrasian equilibrium.
Additionally, the opening prices and closing prices can differ vastly due to
overnight information getting priced into the opening auction; therefore we
opt to remove these overnight returns.
Removing the overnight returns is not a trivial task, and as such this process
is performed as part of the data cleaning process. As a result, algorithm 1, 2
and 3 had to be adjusted to allow for logged-returns as an input, rather than
the original prices 6 . In addition to the overnight return, we also opted to re-
move the first 10 minutes of the continuous trading session i.e. 09:00 to 09:10.
This is because older trading algorithms are still calibrating during this pe-
riod and thus the trades in this period will not be an accurate representation
of the continuous trading session.
5 We highlight the importance of pulling data with the right time zone UTC + 2 for South
African data. We pulled UTC by mistake and thus our trading times are shifted back two
hours.
6 The pair-wise implementation can be found in ftcorr-RealData.R.
30 Chapter 4. Data Engineering

4.3 Creating Data Samples

To investigate the Epps effect, aggregation methods to create various sam-
pling intervals are required. To this end, we will aggregate the TAQ data us-
ing calendar time aggregation methods; specifically the creation of our own
bar data. Additionally, we will perform the novel task of aggregating TAQ
data using intrinsic time aggregation methods; specifically using the frame-
work provided by Derman [14] and our own pragmatic method of aggrega-
tion to highlight the differences between the two estimators. We highlight
for all the aggregation methods, we first obtain the aggregated prices, then
the conversion to returns is performed and finally the overnight returns are
then removed.

4.3.1 Asynchronous Data

The first dataset to create is the cleaned version of the TAQ data, achieved by
aggregating the trades using algorithm 4. Once the aggregation is complete
for each ticker, they get merged into a data frame as shown in figure 4.4 where
the non-trade times are represented with NaNs.

F IGURE 4.4: BTI and NPN Asynchronous Price Sample

The merging is achieved by first pre-populating a data frame with the high-
est available sampling frequency (1 second) over the period of consideration,
then slotting the prices for each asset into their respective times and remov-
ing entire rows of NaNs afterwards.
4.3. Creating Data Samples 31

To create the return matrix required, the returns are computed separately for
each asset over each of the days considered, and the first return for each day
takes on the time index of the second trade in that day, while the first trade
for each day takes on NaN as a placeholder 7 . By computing the returns for
each day separately, we have dealt with the over-night returns. The merging
is then achieved in the same manner as the prices and the result is shown in
figure 4.5.

F IGURE 4.5: BTI and NPN Asynchronous Return Sample

The remaining aggregation methods will be computed using above asyn-

chronous samples 8 .

4.3.2 Calendar Time TAQ Data Aggregation

The calendar time aggregation is simply OHLCV data which is provided by
Bloomberg; however, we opt to create our own bar for the main reason to
append an additional VWAP column as another method to investigate the
Epps effect. The creation of the bar data is done using algorithm 5 9 .
7 afterthe computing the returns
8 This is not an issue for the intrinsic time methods, because the aggregation done here is
VWAP which is the same aggregation used in the intrinsic time samples.
9 The implementation of all the aggregation methods are done using the function gener-

ate_data from TradeDataMain.R.

32 Chapter 4. Data Engineering

Algorithm 5 Bar Data

Require:
1. Ti the unique trading times, i = 1, ..., N
2. Si the aggregated prices, i = 1, ..., N
3. Vi the volume associated with the trade, i = 1, ..., N
4. τ the bar length  . Units of time
Gather trades into their respective bars { J } j , j = 1, ..., M determined by
τ.
for j = 1, ..., M do
if j = 1 then
Set O1 = S1
else
Set O j = Cj−1
end if
Set Hj = max{Si ∈ Jj }
Set L j = min{Si ∈ Jj }
Set Cj = the last Si ∈ Jj
Set Vj∗ = ∑i∈ Jj Vi
Set
∑i∈ Jj Si ∗ Vi
VWAPj =
∑i Vi
Set t∗j = T1 + jτ
end for
return (t∗ = {t∗j } jM=1 , O = {O j } jM=1 , H = { Hj } jM=1 , L = { L j } jM=1 , C =
{Cj } jM=1 , V ∗ = {Vj∗ } jM=1 , VWAP = {VWAPj } jM=1 )

The main difference between the bar data we created and the bar data ex-
tracted from Bloomberg is that Bloomberg’s bar data clocks at exact minutes,
whereas our bar data does not. This difference is due to the fact we wanted
a function that can create bar data for any dataset given. Therefore to avoid
data snooping the first opening price, we began the counter from the time of
the first trade. This was a pragmatic choice because unlike Bloomberg which
will always have a previous closing price to pull from, our finite dataset has
its limitations.
Algorithm 5 10 is presented for 1 asset, however the bar data for the analysis
is for multiple assets. Therefore to make the bar data for more than one asset,
the only point to note is that the T1 used in t∗j = T1 + jτ is computed as the
earliest trade of the day across all the assets across all the days considered 11 .
Then algorithm 5 is applied for each asset.
Figure 4.6 and 4.7 shows the result of algorithm 5 applied to two assets and
converting the closing price and VWAP price to returns respectively.
10 The implementation of algorithm 5 can be found in SynchronousData.R.
11 Forassets which do not have an opening trade the same time as T1 , the opening price is
set to NaN.
4.3. Creating Data Samples 33

F IGURE 4.6: BTI and NPN 10 Minute Closing Bar Return Sam-
ple

F IGURE 4.7: BTI and NPN 10 Minute VWAP Bar Return Sample

A point of detail to note about the creation of figure 4.6 and 4.7 is that the
OHLCV is computed for each asset, and for assets which do not have any
trades within a bar, that row of OHLCV is not computed and therefore skipped.
Returns are then computed for each assets closing and VWAP for each day,
then merged into a data frame in the same manner as the creation of asyn-
chronous returns 12 .
12 We initially made the error for figure 4.7 whereby we computed the OHLCV based on
returns, rather than computing the returns after obtaining the OHLCV for the prices.
34 Chapter 4. Data Engineering

4.3.3 Intrinsic Time TAQ Data Aggregation

Intrinsic time, also known as event time is a method to measure time based
on events rather than the traditional chronological time we humans perceive.
This method of aggregation presents some statistical advantages - namely
the logged-returns are made more gaussian. Additionally, this is the realm
“silicon traders” 13 operate in [29], and they account for a large proportion
of the volume traded [30]. Therefore as a novel application, we will examine
the correlation structure produced from the estimators under this paradigm.

Derman Framework
The first method we will use to aggregate TAQ data in intrinsic time is the
framework provided by Derman [14], where each stock has its own trading
frequency v j . The added benefit from this framework is that it is a natural
way to deal with the asynchrony from high-frequency data through the fact
that each stock has its own trading frequency, and more importantly it pro-
vides an elegant link between intrinsic time and calendar time. However,
this framework has its drawbacks (discussed in section 5.2).
To implement this method of aggregation, we first need the average trades
per day for each stock V̄j over the given data period considered (31/05/2019
- 07/06/2019).

Algorithm 6 Derman Framework

Require:
1. Ti the unique trading times, i = 1, ..., N
2. Si the aggregated prices, i = 1, ..., N
3. Vi the volume associated with the trade, i = 1, ..., N
4. v j = V̄j /(Number of Buckets) the bucket size
Expand the number of observations by repeating each observation Si as
many times as Vi , resulting in I = ∑i Vi observations of Si . Expand such
that the initial ordering of Si is not lost
Set τ = 0
while τv j < I do
τ = τ+1
∀i ∈ [(τ − 1)v j + 1, τv j ], compute

∑i Si ∗ Vi
Pτ =
∑i Vi
end while
return P = { P1 , ..., PM } . M = τ at end of while loop
13 High Frequency Traders
4.3. Creating Data Samples 35

The first point to note about algorithm 6 14 is that it is for one trading day,
and thus the remaining trades at the end of each day which do not have
enough volume to form a bucket are discarded. This choice although devi-
ates from the framework, is justifiable. This is because even though intrinsic
time operates on a separate measurement of time, trading is still performed
on calendar time and at the end of each day the “silicon traders” stop trad-
ing. Furthermore, due to the overnight period, the opening auction can shift
the prices to a completely different level and thus combining the remaining
trades at the end of each day with the first few trades of the next day is not
coherent. Due to this choice, the samples created are not completely syn-
chronous as expected from the framework. This is because V̄j is computed as
the average trades per day over the given period while day to day volume
traded can be different, therefore some assets will have more (less) prices
than the Number of Buckets due to the volume traded in that day being more
(less resp.) than the average. Thus the non-trading times (in intrinsic time)
are filled in with NaNs. Algorithm 6 is computed for each trading day sepa-
rately, then combined afterwards. Finally, the overnight returns are removed
in the same manner as before, by computing the returns for each day sepa-
rately and combining it afterwards. The overnight returns are removed for
the reason that humans operate in calendar time and overnight information
can get priced into the opening auction, therefore changing the price level,
resulting in a return that is not consistent with the continuous trading pro-
cess.
Figure 4.8 shows the resulting intrinsic time return samples for the first trad-
ing day of multiple assets where the average trades per day are computed
over the period (31/05/2019 - 07/06/2019).

F IGURE 4.8: Multiple Ticker Derman Volume Buckets

14 The implementation of algorithm 6 can be found in DermanFrameworkVolumeBuck-
ets.R.
36 Chapter 4. Data Engineering

Lining Up Events
Due to the synchronicity of the Derman framework, the resulting correla-
tion estimates are very similar for both estimators. Thus it is not particularly
meaningful in the context of comparing our estimators. Therefore we employ
another method of aggregating TAQ data in intrinsic time which preserves
the asynchrony of events while retaining the benefit of gaussian returns.

Algorithm 7 Lining Up Events

Require:
j
1. Ti the unique trading times of the jth asset, i = 1, ..., Nj
j
2. Si the aggregated prices of the jth asset, i = 1, ..., Nj
j
3. Vi the volume associated with the trade of the jth asset, i = 1, ..., Nj
4. v = V̄/(Number of Buckets) the bucket size . V̄ is computed from the
most liquid asset
∗
Obtain a sequence of unique trading times { T ∗ }iN∗ =1 across all j assets .
N ∗ ≥ Nj ∀ j
Set τ = 0
for i∗ in T ∗ do
j j j
A j = {Sk }kK=1 . Storage set for the jth asset. Expand Si∗ as many times
j
as Vi∗ and append to the end of the set.
while K~ > ~v do . While any of the storage sets are larger than v
τ = τ+1
for jth asset do
if K j > v then
∀k ∈ [1, v], compute

j ∑ k Sk
Pτ =
v
j
Remove the first v Sk from A j
else
j
Pτ = NaN
end if
end for
end while
end for
j j
return P j = { P1 , ..., PM } . M = τ at end of the loops

The motivation behind algorithm 7 15 is that we want to capture events. This

is because when there is a shock or major event, assets will “synchronise”
in response to the event and once again “de-synchronise” after the event.
Therefore, algorithm 7 works by collecting trades in calendar time until the
bucket is full and a price is printed. Another difference now is that we line up
15 The implementation of algorithm 7 can be found in LiningUpEventsVolumeBuckets.R.
4.3. Creating Data Samples 37

the events based on the most liquid asset which provides a baseline “clock”
for which assets are lined up accordingly.
V̂ in algorithm 7 is computed the same way as V̂j - the average trades per day
over the period considered. The main difference between algorithm 7 is that
the prices need to be computed for all assets at the same time, whereas the
previous algorithms permitted prices for individual assets to be created then
combined. Due to the nature of this aggregation, computing the returns is not
as simple as before, as in we cannot simply use the function diff() anymore.
This is because the function computes x [(1 + lag) : n] − x [1 : (n − lag)], and
therefore if there are no successive prices, the returns will not be computed.
Thus to overcome this issue, we need to first extract the actual prices for
each asset, compute the returns then place them back into their respective
positions.
Removing the overnight returns follows the same process as for algorithm 6,
we apply algorithm 7 for each trading day, then combine them afterwards.
Furthermore, the remaining trades at the end of each day which cannot fill
up a bucket gets discarded - for the same reason we discarded them in the
Derman framework, to focus only on the continuous trading process. How-
ever, this poses an issue for the less liquid stocks when the Number of Buckets
are smaller. This is because the bucket sizes become very large and therefore
some of the less liquid stocks do not have enough trades to form two prices
and thus returns cannot be computed for that day. Due to this issue, we will
ignore the Calendar time equivalent of 1 hour bar samples and focus only on
a bucket frequency of 48 and 480.
Figure 4.9 below illustrates a data sample using algorithm 7 with a bucket
frequency of 48 and the basis ticker is FSR.

F IGURE 4.9: Multiple Ticker Lining Up Events Volume Buckets

From the data samples in figure 4.8 and figure 4.9, we note that there is a
striking difference between these samples. Namely, in figure 4.8, the trades
are near synchronous while in figure 4.9, we have high levels of asynchrony
- allowing us to gain further insights into the two estimators and how they
compare.
 39

Chapter 5

Data Science

Turning our attention to real financial data, we perform the novel application
by studying the Epps effect through various methods of aggregating Trade
and Quote (TAQ) data. Specifically, we compare calendar time based sam-
pling with volume time sampling methods.

5.1 Calendar Time

For the calendar time sampling methods we focus on the closing prices and
Volume Weighted Average Price for different sized bar data.

5.1.1 Closing Prices

Using the TAQ data, we create 1 minute, 10 minute and 1 hour OHLCV bar
data for the period. From the Monte Carlo experiment, we know that the two
estimators produce the same results when trades are synchronous; however,
some of the less liquid stocks do not have any trades within a given bar and
therefore the bar data created is not truly synchronous. Thus the MM and HY
estimates will differ slightly with HY having a higher estimate. We apply the
MM and HY estimators to the closing prices of the various bar data along
with the raw TAQ data to see the effect of increasing the sampling frequency
to its highest available frequency.
Figure 5.1 1 2 (a) though to (d), the MM estimator is applied to 1 hour, 10
minute, 1 minute closing bar data and TAQ data respectively. Figure 5.1 (e)
through to (h), the HY estimator is applied to 1 hour, 10 minute, 1 minute
closing bar data and TAQ data respectively. On each of the plots, the aver-
age magnitude of the correlation is shown. It is clear that the Epps effect is
present for both estimators when considering the closing price samples. This
shows that the Epps effect is present even at much higher sampling frequen-
cies that what [4] originally considered.
1 Figure 5.1 can be reproduced using Closing.R and TAQ.R.
2 Figure 5.2 can be reproduced using VWAP.R and TAQ.R.
40 Chapter 5. Data Science

(a) MM 1 HR Close Price (e) HY 1 HR Close Price

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.2829 −0.5 AGL |ρij| = 0.2912 −0.5

NPN NPN
BTI BTI
−1.0 −1.0

MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI
(b) MM 10 Min Close Price (f) HY 10 Min Close Price
ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.2334 −0.5 AGL |ρij| = 0.2347 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

(c) MM 1 Min Close Price (g) HY 1 Min Close Price

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.1384 −0.5 AGL |ρij| = 0.1466 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

(d) MM TAQ (h) HY TAQ

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.0355 −0.5 AGL |ρij| = 0.1051 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

F IGURE 5.1: Investigating the Epps effect by aggregating TAQ

data into closing bar prices. From (a) to (d), we have the MM
estimator applied to 1 hour, 10 minute, 1 minute closing bar
data and TAQ data respectively using algorithm 2. From (e) to
(h), we have the HY estimator applied to 1 hour, 10 minute, 1
minute closing bar data and TAQ data respectively using algo-
rithm 3. The Epps effect [4] is demonstrated with both estima-
tors and it persists with the MM but only slightly for the HY
estimator on the TAQ data.
5.1. Calendar Time 41

What is more interesting is that when the sampling frequency increases to

the highest available frequency, the correlation from the MM estimates drops
towards zero as expected from the Epps effect, but for the HY estimates, the
correlation from the TAQ data does not die out completely and stays near
the correlation estimates from the 1 minute bars. Indicating that the HY es-
timator is manufacturing correlation and thereby not well suited for high-
frequency data due to the upward bias. This is because the Epps effect is
very real and present in both estimators.
Additionally, the drop in correlation from figure 5.1 (e) through to (g) goes
against the claim of Hayashi and Yoshida that the Epps effect is a bias in the
estimator for which their estimator is immune to, as we have clearly demon-
strated otherwise. Even though their estimator does seem to be immune to
the Epps effect arising from asynchrony, it does not seem to always be the
case as seen in figure 3.5 (f), therefore given that there is no seamless way to
decompose the various factors contributing towards the Epps effect it is hard
to pinpoint what effect the multiple contributions are having with regards to
the Epps effect.
The correlation estimates make quite a lot of sense, the top right corner for the
various sub-figures in figure 5.1 are highly correlated, and these tickers are all
from the banking sector, therefore validating the fact that we have performed
the analysis correctly and algorithm 2 and 3 are implemented correctly.
Another interesting note is as the sampling intervals become smaller, the
correlations drop, but become positively correlated. This is clearly demon-
strated in figure 5.1 (a) to (b) and (e) to (f). Suggesting the argument of asyn-
chrony and lead-lags may not be sufficient in fully explaining the entirety of
the Epps effect, but rather there is an underlying change in correlation struc-
ture depending on the sampling interval.
Using closing prices for the various bar data discards a lot of the information
given by the financial market as it is just another sample from the available
sample of TAQ data. Thus we turn our attention to another method to create
bar data, which encapsulates more of the information given within a bar.

5.1.2 Volume Weighted Prices

Looking at another method of aggregating TAQ data, we look towards the
Volume Weighted Average Prices (VWAP). This representation of bar data
is better than the closing prices because the VWAP includes the information
from all trades within a given bar by means of averaging. Using this aggrega-
tion method, we apply the MM and HY estimators to the 1 hour, 10 minute,
1 minute and TAQ data to see how this compares against the closing prices.
From figure 5.2 (a) to (d), the MM estimator is applied to the 1 hour, 10
minute, 1 minute VWAP bar data and TAQ data respectively; and from figure
5.2 (e) to (h), the HY estimator is applied to the 1 hour, 10 minute, 1 minute
VWAP bar data and TAQ respectively. Comparing the two estimators, the
42 Chapter 5. Data Science

Epps effect persists with the MM estimator into the TAQ data, while the cor-
relation does not completely die out with the HY estimator as we go further
into the high-frequency spectrum - indicating that the HY is manufacturing
correlations through its multiple contributions. Further attesting to our ar-
gument that the HY estimate is biased for high-frequency event data.
The correlation structures in figure 5.2 is very similar to that of 5.1, the bank-
ing sectors are still strongly positively correlated. The main difference be-
tween the two correlation structures is that the VWAP aggregation seems to
accentuate the existing correlation structure in the Closing aggregation. This
is very interesting because the aggregation methods are quite different, the
VWAP incorporates more information within the given bar by means of aver-
aging while the closing prices are mere samples from the finite TAQ sample,
yet their correlation structures remained very similar. Additionally, the Epps
effect is once again present in this method of aggregation - once again sug-
gesting there seems to be a structural change in correlation that is dependent
on the sampling interval.
An interesting point to note in figures C.3 and C.4, is that we initially made
the error of computing the OHLCV based on returns, rather than computing
the returns after obtaining the OHLCV prices. Therefore, initially we had
the closing and VWAP returns computed on the TAQ data for the various
intervals, instead of computing the returns over the various intervals. The
interesting thing about this is that we still saw an Epps effect under these
circumstances, where the returns are computed from the highest available
sampling frequency, rather than dependent on the various sampling inter-
vals. Therefore the Epps effect was also inadvertently achieved by sampling
the TAQ returns at various sampling frequencies. This begs the question as
to what the Epps effect truly is. How and why did it still show up even when
returns are computed at the highest available sampling frequency?
The VWAP bars share the same issue as the closing bars; some of the less
liquid stocks will not have any trades within a given bar, thus the bar data is
not truly synchronous and the two estimators differ slightly. It must be noted
that although the VWAP incorporates more information from the stocks it
has its own issues; namely that it hides any jumps the price paths may have
into the average, therefore acting as a smoothing operator similar to that
of a Moving Average. An additional issue is that the aggregation between
the bars lacks consistency, this is because aggregating the data in calendar
time means that different bars will be averaged with different volume sizes.
Therefore a more suitable way to incorporate information from the price path
is to look at intrinsic time aggregation.
5.1. Calendar Time 43

(a) MM 1 HR VWAP (e) HY 1 HR VWAP

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.3172 −0.5 AGL |ρij| = 0.3201 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI
(b) MM 10 Min VWAP (f) HY 10 Min VWAP
ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.256 −0.5 AGL |ρij| = 0.2441 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

(c) MM 1 Min VWAP (g) HY 1 Min VWAP

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.1269 −0.5 AGL |ρij| = 0.145 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

(d) MM TAQ (h) HY TAQ

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.0355 −0.5 AGL |ρij| = 0.1051 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

F IGURE 5.2: Investigating the Epps effect by aggregating TAQ

data into VWAP bar prices. From (a) to (d), we have the MM
estimator applied to 1 hour, 10 minute, 1 minute closing bar
data and TAQ data respectively using algorithm 2. From (e) to
(h), we have the HY estimator applied to 1 hour, 10 minute, 1
minute closing bar data and TAQ data respectively using algo-
rithm 3. The Epps effect [4] is demonstrated with both estima-
tors and it persists with the MM but only slightly for the HY
estimator on the TAQ data.
44 Chapter 5. Data Science

5.2 Intrinsic Time

Intrinsic time is a different paradigm where time is now measured in terms
of machine time - the clock ticks based on events rather than what humans
perceive as time. Working in volume time presents significant statistical ad-
vantages. The benefit is that it allows the partial recovery of Normality and
the IID assumption. Additionally, sampling in a volume clock metric deals
with the issue of random and asynchronous trade data [29]. Furthermore,
under this paradigm, there is a framework presented by Derman [14] which
provides an intuitive method to synchronise the TAQ data across the various
assets.

5.2.1 Derman Framework

Derman assumed that each stock has their its intrinsic time scale which is
constant through time. He defines the stocks’ trading frequency v j as the
number of intrinsic time ticks that occur for one calendar second [14].
He gives the relationship between the flow of calendar time t and the flow of
intrinsic time τj as

dτj = v j dt. (5.2.1)

More importantly, he was able to show that the correlation in intrinsic time
πij is the same as the correlation in calendar time ρij [14]. This framework
provides a method to create synchronous price paths in intrinsic time which
will allow the recovery of the correlation in calendar time. However, there
is currently no framework that provides a method to create price paths in
intrinsic time while allowing for the sampling intervals (in intrinsic time) to
reduce down to each individual volumes. Thus we are unable to study the
equivalent of TAQ data in intrinsic time.
To study the Epps effect with this framework, we have to alter Derman’s
methodology slightly. Instead of assuming v j as the number of trades per cal-
endar second, we assume v j to be the number of trades per unit of sampling
interval considered. Although we altered his method slightly, the maths
showing that the correlations are dimensionless and independent of the var-
ious time measurements still holds.
Therefore to apply this sampling scheme, bucket sizes must be chosen for
each stock. This is determined by the rough equivalent bar length in calendar
time (i.e. to create the equivalent of 10 min calendar time bars in intrinsic
time, we divide the average volume per day by 48). Since the bucket sizes
are computed from the average volume per day, the price paths will not be
fully synchronous due to different volume amounts traded per day. Some
stocks will have more (less) volume buckets if the trades for the day are above
(below) the average for the stock (respectively).
5.2. Intrinsic Time 45

(a) MM 8 Volume Bucket (d) HY 8 Volume Bucket

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.2209 −0.5 AGL |ρij| = 0.266 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI
(b) MM 48 Volume Bucket (e) HY 48 Volume Bucket
ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.0879 −0.5 AGL |ρij| = 0.0914 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

(c) MM 480 Volume Bucket (f) HY 480 Volume Bucket

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.023 −0.5 AGL |ρij| = 0.0338 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

F IGURE 5.3: Using the Derman framework [14] to compare the

two estimators. From (a) to (c), we have the MM estimator ap-
plied to the calendar time equivalent of 1 hour, 10 minute and 1
minute bar data respectively using algorithm 2. From (d) to (f),
we have the HY estimator applied to the calendar time equiv-
alent of 1 hour, 10 minute and 1 minute bar data respectively
using algorithm 3. The Epps effect is clearly present for both
estimators under the paradigm of intrinsic time.

This highlights the first issue with the framework. The assumption that each
stocks’ intrinsic time scale is constant through time is a strong assumption
to make, where the validity is questionable. This is because the trading fre-
quency changes over time depending on factors such as time of day or rele-
vant news reports. For instance, if news comes out that a company is about
to undergo liquidation, traders will try to square their positions therefore in-
creasing the trading frequency. The second issue is that (5.2.1) is the key to
linking the correlation in intrinsic time to the correlation in calendar time, but
46 Chapter 5. Data Science

this assumption induces a continuity assumption from the calendar time into
the intrinsic time given the linear relationship which re-scales the time given
by standard continuous-time stochastic processes. It must therefore be noted
that it does not fully achieve the effect of converting stocks into intrinsic time
- where time ticks purely on the events. Elaborating further on why this does
not achieve the full conversion into intrinsic time is because a unit interval
in intrinsic time can be thought of as a stochastic interval in calendar time,
where the stopping rule is determined by the number of trades counted.
Using parts of algorithm A.1. from [30] in algorithm 6, we create the 1 hour,
10 minute and 1 minute calendar time equivalent data by using 8, 48 and 480
buckets per day respectively to create the intrinsic time samples.
From figure 5.3 3 (a) through to (c), we have the MM estimator applied to
the calendar time equivalent of 1 hour, 10 minute and 1 minute bar data
respectively using algorithm 2. From figure 5.3 (d) through to (f), we have
the HY estimator applied to the calendar time equivalent of 1 hour, 10 minute
and 1 minute bar data respectively using algorithm 3. It is clear that the Epps
effect still exists under this completely different method of aggregating TAQ
data.
This is interesting because this shows that the Epps effect does not only ex-
ist in the paradigm of calendar time, it is also present under the event time
paradigm. What is even more interesting is that the correlation structures
change depending on the sampling interval used, indicating that correlations
are not indeed dimensionless as suggested by Derman [14], and that the Epps
effect seems to be intrinsically linked to the sampling intervals chosen, there-
fore further attesting to the idea that the Epps effect cannot be fully explained
with asynchrony or lead-lags. In addition, the correlation structure in figure
5.3 is very different to that of figure 5.1 and 5.2, indicating that the correla-
tions are not preserved across these various measurements of time as Derman
proved. Finally, this method faces similar issues to that of the VWAP where
the jumps are hidden into the averages. Although we have highlighted a few
of the pitfalls regarding this framework, this is still the most seamless frame-
work provided in the literature which ties together the ideas from intrinsic
time to calendar time.
This however does not answer the main question as to which estimator is
the more efficient of the two; simply because the aggregation method creates
data which is very close to being synchronous, therefore the two estimators
behave extremely similarly as seen in figure 5.3. To this end, by employing
the ideas from the intrinsic time framework [14], we create our own method
of aggregating TAQ data, specifically focused on determining how the two
estimators differ.
3 Figure 5.3 can be reproduced using Derman.R.
5.2. Intrinsic Time 47

5.2.2 Lining Up Events

For this method of aggregation, the averaging method used in the Derman
framework is also used here. However the difference is that we do not as-
sume that each stock has its own intrinsic time. What we do is we count the
events in terms of calendar time, and each stock “collects” trades in calendar
time until the bucket is full and a price is printed 4 . The bucket sizes are deter-
mined by the average trades per day of the most liquid stock which provides
the baseline “clock” for which all the assets are lined up accordingly.
The rationale behind this method of aggregation is that we want the prices
to line up when there is an event triggering more trades. This is because
when there is a large event, trades across various assets will “synchronise”
in response to the event and once again “de-synchronise” after the event has
passed. There is no framework linking this method of aggregation to the
correlation in calendar time; however the advantages this method possesses
is that by employing the volume time aggregation method we benefit from
have log-returns that are more gaussian, in addition, this is a natural way
to aggregate events in an asynchronous manner. Finally, since this method
does not resolve the issue of asynchronous trading times, in the context for
comparing the two estimators, this is exactly what we want - a clear example
to demonstrate how the two estimators differ.
An issue was encountered when creating the 1 hour equivalent calendar time
sample. Due to the large bucket sizes based on the most liquid stock, the
stocks with low liquidity simply did not have enough trades in a given day
to fill up more than two buckets. Since there is no clear method to fix the de-
fective sample, we discard this sample and focus only on the 10 minute and
1 minute calendar time equivalent samples. This does not affect the analysis
since correlations on short time scales is what is of interest.
Figure 5.4 5 (a) and (b) is the MM estimator applied to the 10 minute and
1 minute calendar time equivalent with algorithm 2 and (c) and (d) is the
HY estimator applied to the 10 minute and 1 minute calendar time equiva-
lent with algorithm 3. The first thing to notice is that even under this asyn-
chronous sampling scheme the Epps effect is still present with both estima-
tors. The second thing to notice is how different the two estimators are: the
MM has correlation near zero while HY has an extremely high correlation.
The 1 minute equivalent calendar bar (d) has a significantly higher correla-
tion when compared to the other aggregation methods - indicating that HY
is manufacturing correlation.
Figure 5.4 is the key figure in this paper that demonstrates the HY is bi-
ased for computing the correlation of events. This is because HY does not
treat missing observations as proper missing observations, it compensates
the missing observations through the multiple contributions to bring up the
estimate and therefore manufacturing correlation. The MM treats missing
4 The time steps used here is still intrinsic time
5 Figure 5.3 can be reproduced using LiningUpEvents.R.
48 Chapter 5. Data Science

observations as truly missing observations, it does not compensate the cor-

relation but rather uses a lossless interpolation between events, therefore
it does not overestimate the correlation. Therefore further showcasing our
preference of the MM estimator over the HY estimator for high-frequency
trading - where the events are truly discontinuous and asynchronous. The
HY estimator is most certainly the better estimator of the two, given that
the observed prices in the financial market are samples from an underlying
continuous-time stochastic process; however the MM estimator is the better
estimator of the two, given that the high-frequency finance world is truly
discontinuous, discrete and asynchronous events.

(a) MM 48 (c) HY 48
ρ ρ
FSR FSR
1.0 1.0
SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0
SOL SOL 0.0

MNP MNP
AGL |ρij| = 0.0073 −0.5 AGL |ρij| = 0.4265 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

(b) MM 480 (d) HY 480

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.0046 −0.5 AGL |ρij| = 0.2996 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

F IGURE 5.4: Comparing the two estimators using a sampling

method to line up events. (a) and (b) is the MM estimator ap-
plied to the 10 minute and 1 minute calendar time equivalent
with algorithm 2, (c) and (d) is the HY estimator applied to the
10 minute and 1 minute calendar time equivalent with algo-
rithm 3. The Epps effect is once again present, but we see that
HY is manufacturing correlation.
 49

Chapter 6

Concluding Remarks

In this report, we have demonstrated that the MM and HY estimators differ

under asynchronous conditions which is the exact conditions we are inter-
ested in addressing with regards to high-frequency finance. Therefore al-
though we have argued for the efficacy of the MM estimator over the HY
estimator from a data-informed view, this report does not concretely show
which estimator is the better of the two. The whole argument as to why the
MM estimator is more believable lies in the fact that it possesses lossless in-
terpolation between event data and it correctly captures the existence of the
Epps effect and the fact that it is currently unclear what the multiple contri-
butions are truly doing with regards to the Epps effect in the HY estimator.
The reason why one cannot definitively show that the MM estimator per-
forms better is because the two estimators only differ under tick-by-tick asyn-
chronous trade data which is also the exact case where we want to deter-
mine the more effective estimator. Furthermore, given that there is currently
no cohesive framework that ties together the data generating process and
the continuous-time stochastic process; no seamless framework that ties up
all the loose ends regarding the Epps effect; no satisfactory method to syn-
chronise the TAQ data without interpolation and no method to falsify high-
frequency finance being samples from a continuous-time stochastic process -
one is left with very limited methods to definitively demonstrate which es-
timator performs better. What we have demonstrated in this report is that
there is some inconsistency in the HY estimator when it comes to investigat-
ing the Epps effect and that the estimator seems to be manufacturing corre-
lation which does not exist.
This report does not show anything new that was previously unknown about
the Epps effect. What this report does show is the existence of the Epps ef-
fect on the JSE through various methods of aggregating TAQ data, which
suggests that the Epps effect explained purely by asynchrony and lead-lag is
insufficient, further validating Thomas Epps’ initial thought that there seems
to be an underlying change in the correlation structure which is dependent
on the sampling intervals used [4].
A natural extension to what we have done here is to increase the computation
speed of algorithm 2 using Fast Fourier techniques. Specifically, design a Fast
Fourier Transform (FFT) version of algorithm 2 for the synchronous case [31]
50 Chapter 6. Concluding Remarks

and a Non-Uniform FFT (NUFFT) version for the more general asynchronous
case.
Additional extensions regarding the Monte Carlo experiments in section 3
include the multivariate Hawkes process [25] method of generating asyn-
chronous data and to study the effect of the correlation estimates by chang-
ing the various parameters pertaining to the Hawkes process calibration and
simulation. Another extension on the Monte Carlo side will be to recover the
results from [6] and firstly see if the MM estimator recovers the same esti-
mates as their Fourier method, and secondly use these results to decompose
the Epps effect into the component arising from asynchrony and the remain-
ing factors.
Finally, all the code listing used in this report can be found on Github [32] and
the steps outlining the recovery of the results can be found in the [Link]
document to ensure straightforward replication of the results.
 51

Appendix A

Supporting Algorithms

Algorithm 8 Re-scale Trading Times

Require:
1. (n x m) matrix T of asynchronous sampled times
tmin ← minimum value of T
tmax ← maximum value of T
for j = 1 to m do
for i = 1 to n do
2π (t −t )
τij ← tmaxj −t min
min
end for
end for
return (τ)

Algorithm 9 Kanatani Weight

Require:
1. (n x 2) matrix τ and P of re-scaled asynchronous sampled times and
prices
I. Initialize W matrix
W ← 0 [Ni x Nj ] matrix of 0’s
for j = 1 to Ni do
for i = 1 to Nj do
j j
if (tik−1 , tik ] ∩ (tl −1 , tl ] 6= ∅ then
wkl ← 1
end if
end for
end for
return (W)

12

1 Algorithm 8 can be found in ftcorr.R as an auxiliary function and was provided by [18].
2 Algorithm 9 can be found in ftcorr.R as an auxiliary function.
52 Appendix A. Supporting Algorithms

Algorithm 10 Simulating Geometric Brownian motion

Require:
1. n number of price points to simulate
2. µ (d x 1) vector of drift parameters
3. Σ (d x d) covariance matrix
4. start price (d x 1) vector of S(0)
Procedure for the ith asset:
1. generate Z ∼ Nd (0, Idxd )
√
2. set Si (tk+1 ) = Si (tk ) exp (µi − 21 σi2 )(tk+1 − tk ) + tk+1 − tk ∑dk=1 Aik Zk


return (S)

Subject to the condition S(0) = start price and A is the Cholesky decomposi-
tion of Σ. Algorithm 10 is provided by [22] 3 .

Algorithm 11 Simulating Merton Model

Require:
1. n number of price points to simulate
2. µ (d x 1) vector of drift parameters
3. Σ (d x d) covariance matrix
4. λ (d x 1) vector of the Poisson process parameter
5. a (d x 1) vector of lognormal location parameter
6. b (d x 1) vector of lognormal standard deviation
7. start price (d x 1) vector of S(0)
Procedure for the ith asset:
1. generate Z ∼ Nd (0, Idxd )
2. generate Ni ∼ Poisson(λi (tk+1 − tk ))
3. generate Z2 ∼ N1√ (0, 1)
4. set M = ai Ni + bi Ni Z2
5. set Xi ( t k + 1 ) = Xi (tk ) + (µi − 12 σi2 )(tk+1 − tk ) +
√
tk+1 − tk ∑dk=1 Aik Zk + M
6. S = exp( X )
return (S)

Subject to the condition X (0) = ln(start price) and A is the Cholesky decom-
position of Σ. Algorithm 11 is provided by [22] 4 .
3 Algorithm 10 can be found in GBM.R.
4 Algorithm 11 can be found in Merton Model.R.
Appendix A. Supporting Algorithms 53

Algorithm 12 Simulating GARCH(1,1)

Require:
1. n number of price points to simulate
2. θ (2 x 1) mean reverting rate
3. λ (2 x 1)
4. w (2 x 1) vector long term variance
5. ρ correlation
6. starting variance (2 x 1) vector of starting variance
7. starting price (2 x 1) vector of starting price
Procedure for the ith asset:
1. generate Z ∼ N (0, 1)
2. set σi2 (tk+1 ) = σi2 (tk ) + θi (w1 − σi2 (tk ))(tk+1 − tk ) +
p
2λi θi (tk+1 − tk )σi (tk ) Z
3. create Σ∗ correlation matrix based on σ2 (tk+1 )
4. generate Z ∗ ∼ Nd (0, Idxd )
√
5. set Xi (tk+1 ) = Xi (tk ) + tk+1 − tk ∑dk=1 Aik ∗ Z∗
k
6. S = exp( X )
return (S)

Subject to the condition X (0) = ln(start price), σ(0) = starting variance and
A∗ is the Cholesky decomposition of Σ∗ 5 .

Algorithm 13 Simulating Variance Gamma

Require:
1. n number of price points to simulate
2. µ (d x 1) vector of drift parameters
3. Σ (d x d) covariance matrix
4. β (d x 1) scale parameter of Gamma
5. start price (d x 1) vector of S(0)
Procedure for the ith asset:
1. generate Yi ∼ Gamma(tk+1 − tk /β i , β i )
2. generate Z ∼ Nd (0, Idxd )
√
3. set Xi (tk+1 ) = Xi (tk ) + µYi + Yi ∑dk=1 Aik Zk
return (X)

Subject to the condition X (0) = start price and A is the Cholesky decompo-
sition of Σ. Algorithm 13 is provided by [22] 6 .
5 Algorithm 12 can be found in GarchAndersen.R, while the specification from [12] can be

found in GarchReno.R.
6 Algorithm 13 can be found in Variance Gamma.R.
54 Appendix A. Supporting Algorithms

Algorithm 14 Simulating Ornstein Uhlenbeck

Require:
1. n number of price points to simulate
2. µ (d x 1) vector of long term prices
3. Σ (d x d) covariance matrix
4. θ (d x 1) vector of mean reverting rate
5. start price (d x 1) vector of starting prices
Procedure for the ith asset:
1. generate Z ∼ Nd (0, Idxd )
2. set Xi (tk+1 ) = Xi (tk ) + θi (ln(µi ) − Xi (tk ))(tk+1 − tk ) +
√ d
tk+1 − tk ∑k=1 Aik Zk
3. S = exp( X )
return (S)

Subject to the condition X (0) = ln(start price) and A is the Cholesky decom-
position of Σ 7 .

7 Algorithm 14 can be found in Ornstein Uhlenbeck.R.

 55

Appendix B

Appendix Derivation

B.1 Proof for Theorem 2.1.1 and 2.1.2

B.1.1 Proof for aq (Σ)
The proof is based on [1] and [15], and we complete the additional results
that were omitted in their papers.
Assume the Bachelier paradigm and further suppose the price process fol-
lows a one-dimensional Itô process

dp(t) = β(t)dt + σ (t)dW. (B.1.1)

Remark B.1.1 The price process need not be a one-dimensional Itô process, we as-
sume so for ease of derivation. The price process in general follows (2.1.1)
Furthermore, we can assume the drift term has no contribution and can there-
fore ignore it. Thus the price process becomes

dp(t) = σ(t)dW. (B.1.2)

This is acceptable because Malliavin and Mancino show that the contribution
for the drift term is zero [1]. Additionally Malherbe argues that ignoring the
drift term implies an efficient market [15].
We now introduce the Gaussian variables

Gk := ak (dp), Gk0 := bk (dp).

From Corollary (3.2.6) to the Martingale Representation theorem in [33], we

know that Gk and Gk0 is a martingale. Therefore, Gk and Gk0 are Gaussian
variables with E( Gk ) = E( Gk0 ) = 0 [15].
We now calculate the covariance of the Gaussian variables
56 Appendix B. Appendix Derivation

E( Gk Gl ) = E( ak (dp) al (dp))
h 1 Z 2π Z 2π i (B.1.3)
=E 2 cos(kt)σ(t)dW cos(lt)σ (t)dW .
π 0 0

Now using the Itô isometry and (2.1.2), (B.1.3) becomes

Z 2π
1
E( Gk Gl ) = 2 Σ(t) cos(kt) cos(lt)dt. (B.1.4)
π 0

Using the identity

1
cos(kt) cos(lt) = (cos(k − l )t + cos(k + l )t),
2
(B.1.4) becomes

1
Z 2π
E( Gk Gl ) = Σ(t) cos(kt) cos(lt)dt
π2 0
Z 2π
1
= Σ(t)[(cos(k − l )t + cos(k + l )t)]dt
2π 2 0
Z 2π Z 2π
1 1
= Σ ( t )( cos ( k − l ) tdt + Σ(t)(cos(k + l )tdt
2π 2 0 2π 2 0
1
= (a (Σ) + ak+l (Σ)).
2π |k−l |
(B.1.5)
The energy identity 1 is

||Σ||2L2 = ∑( ak (Σ))2 + (bk (Σ))2 . (B.1.6)

k

Furthermore,

∑(ak (Σ))2 ≤ ||Σ||2L2 . (B.1.7)

k

q
Now for q > 0 consider the random variable UN as the discrete convolution
of the Gaussian variables, where

N
1
∑ Gk Gk+q .
q
UN := (B.1.8)
N k =1

Using (B.1.5), we get that

1 Better known as Parseval’s theorem
B.1. Proof for Theorem 2.1.1 and 2.1.2 57

N
1
∑ E(Gk Gk+q )
q
E (U N ) =
N k =1
N
1 1
=
N ∑ 2π (a|k−k−q| (Σ) + ak+k+q (Σ))
k =1
N
1 1
=
N ∑ 2π (aq (Σ) + a2k+q (Σ)) (B.1.9)
k =1
N
1  1
aq (Σ) + ∑ a2k+q (Σ)

=
2π N k =1
N
1
∑ a2k+q (Σ)
q
2πE(UN ) = aq (Σ) +
N k =1
= aq (Σ) + R N .

We note that Malliavin and Mancino omitted the 2π which is necessary in

order to recover (2.1.1).
q
We now show that 2πE(UN ) → aq (Σ) as N → ∞ by using the Cauchy-
Schwarz inequality

1 N 1 N 2 1 N 2 1
Rn = | ∑ a2k+q (Σ)| ≤ ( ∑ 1 ) 2 ( ∑ a2k+q (Σ)) 2
N k =1 N k =1 k =1
N
1 1
= √ ( ∑ a22k+q (Σ)) 2 (B.1.10)
N k =1
1
≤ √ ||Σ|| L2 .
N

Therefore, R N → 0 as N → ∞ 2 .
q q
We now want to show that lim N →∞ UN = E(UN ) in L2 . We first compute

1
∑0
q
E[(UN )2 ] = E( Gk2 Gk20 +q ). (B.1.11)
N2 0≤k,k ≤ N

and consider an R2 -valued normal variable ( G1 , G2 ) and denote

λi := E( Gi2 ), µ := E( G1 G2 ),

and define

µ
Z := G2 − G1 .
λ1
2
∑kN=1 a22k+q (Σ)) ≤ ∑k a2k (Σ)
58 Appendix B. Appendix Derivation

We note that E( ZG1 ) = E( G1 G2 − λ1 G12 ) = E( G1 G2 ) − λ1 E( G12 ) = E( G1 G2 ) −

µ µ

E( G1 G2 ) = 0 and therefore G1 and Z are independent. Additionally, since

E( G1 ) = E( G2 ) = 0, we get 3

µ µ2
E( G12 G22 ) = E[ G12 ( Z2 + 2 G1 Z + 2 G12 )]
λ1 λ1
µ µ2
= E( G12 Z2 ) + 2 E( G13 Z ) + 2 E( G14 )
λ1 λ1
µ2
= E( G12 Z2 ) + E( G14 )
λ21
µ2
= E( G12 ) E( Z2 ) + E( G14 )
λ21
µ µ2 µ2
= E( G12 ) E( G22 − 2 G1 G2 + 2 E( G12 )) + 2 E( G14 )
λ1 λ1 λ1
µ µ2 µ2
= E( G12 ) E( G22 ) − 2 E( G12 ) E( G1 G2 ) + 2 E( G12 ) E( G12 ) + 2 E( G14 ).
λ1 λ1 λ1
(B.1.12)
Furthermore, since Gk is a Gaussian variable with mean 0. We know the
fourth uncentered moment is E( Gk4 ) = 3E( Gk2 )2 . Therefore, (B.1.12) becomes

µ µ2 µ2
E( G12 G22 ) = E( G12 ) E( G22 ) − 2 E( G12 ) E( G1 G2 ) + 2 E( G12 ) E( G12 ) + 2 E( G14 )
λ1 λ1 λ1
µ µ2
= E( G12 ) E( G22 ) − 2 E( G12 ) E( G1 G2 ) + 4 2 E( G12 ) E( G12 )
λ1 λ1
µ µ2
= E( G12 ) E( G22 ) − 2 λ1 µ + 4 2 λ1 λ1
λ1 λ1
= E( G12 ) E( G22 ) + 2µ2
= E( G12 ) E( G22 ) + 2E( G1 G2 )2 .
(B.1.13)
Putting this together, we get
3 E( G3 Z ) = 0 due to independence and that the third uncentered moment E( G13 ) = 0
1
B.1. Proof for Theorem 2.1.1 and 2.1.2 59

q q
E[(UN − E(UN ))2 ]
q q q q
= E[(UN )2 − 2UN E(UN ) + ( E(UN ))2 ]
q q
= E((UN )2 ) − [ E(UN )]2
1 1
= 2 ∑ E( Gk2 Gk20 +q ) − a q ( Σ )2
N 0≤k,k0 ≤ N (2π )2
1 1
=
N2 ∑0 [ E(2k ) E( Gk20 +q ) + 2( E( Gk Gk0 +q ))2 ] −
(2π )2
a q ( Σ )2
0≤k,k ≤ N
1 1 1
=
N2 ∑0 [ E(2k ) E( Gk20 +q ) +
2π 2
( a|k−k0 +q| (Σ) + ak+k0 +q (Σ))2 ] −
(2π )2
a q ( Σ )2
0≤k,k ≤ N
1 1 1
∑0
q
= ( E(UN ))2 + [ ( a|k−k0 +q| (Σ) + ak+k0 +q (Σ))2 ] − a q ( Σ )2
N2 0≤k,k ≤ N
2π 2 (2π )2
1
2 1 1 1
=
2π
aq (Σ) + 2
N ∑0 [
2π 2
( a|k−k0 +q| (Σ) + ak+k0 +q (Σ))2 ] −
(2π )2
a q ( Σ )2
0≤k,k ≤ N
1
=
2π 2 N 2 ∑0 ( a|k−k0 +q| (Σ) + ak+k0 +q (Σ))2
0≤k,k ≤ N
1
≤ ||Σ||2L2
N
(B.1.14)
As N → ∞, N || Σ || L2
1 2 → 0. Thus it follows that

q q
lim UN = E(UN ) in L2 .
N →∞

q
Coupled with the fact that 2πE(UN ) → aq (Σ) as N → ∞, we get

q q
2πUN → 2πE(UN ) → aq (Σ),

in probability as a natural consequence from L p convergence as N → ∞.

Therefore,

N
2π
∑ (as (dp)as+q (dp)),
q
aq (Σ) = lim 2πUN = lim ∀q > 0.
N →∞ N →∞ N
s =1

aq (Σ) has been proved. The remaining univariate and multivariate cases will
be less rigorous and more focused on achieving the correct scaling factors.

B.1.2 Proof for a0 (Σ)

We consider E(bs2 (dp)). Using the Itô isometry we get,
60 Appendix B. Appendix Derivation

h 1 Z 2π Z 2π i
E(bs2 (dp)) =E 2 sin(st)σ(t)dW sin(st)σ(t)dW
π 0 0
Z 2π (B.1.15)
1
= 2 sin2 (st)Σ(t)dt.
π 0

Using the identity sin2 (st) = ( 12 − 12 cos(2st)), (B.1.15) becomes

1 2π
Z
E(bs2 (dp)) = 2 sin2 (st)Σ(t)dt
π 0
Z 2π Z 2π (B.1.16)
1 1
= Σ ( t ) dt − cos(2st)Σ(t)dt.
2π 2 0 2π 2 0

Using the identity cos(2st) = 2 cos2 (st) − 1, (B.1.16) becomes

1 2π
Z
1 2π Z
E(bs2 (dp)) = Σ(t)dt − cos(2st)Σ(t)dt
2π 2 0 2π 2 0
Z 2π Z 2π Z 2π
1 1 1
= Σ(t)dt − cos (st)Σ(t)dt +
2
Σ(t)dt.
2π 2 0 π2 0 2π 2 0
(B.1.17)
2π
Using the fact that E( a2s (dp)) = π12 0 cos2 (st)Σ(t)dt, we get that
R

2
a0 (Σ) = E( a2s (dp)) + E(bs2 (dp)). (B.1.18)
π
1
Thus we get a scaling factor of 2 which is not present in [1].

B.1.3 Proof for bq (Σ)

We consider E( as (dp)bs+q (dp)). Using the Itô isometry we get,

h 1 Z 2π Z 2π i
E( as (dp)bs+q (dp)) = E cos ( st ) σ ( t ) dW sin ( s + q ) tσ ( t ) dW
π2 0 0
Z 2π
1
= 2 cos(st) sin(s + q)tΣ(t)dt.
π 0
(B.1.19)
1
Using the identity sin(lt) cos(st) = 2 [sin( l + s)t + sin(l − s)t, (B.1.19) be-
comes
B.1. Proof for Theorem 2.1.1 and 2.1.2 61

1 2π Z
1
Z 2π
E( as (dp)bs+q (dp)) = 2
sin(qt)Σ(t)dt + sin(2s + q)tΣ(t)dt
2π 0 2π 2 0
1
= (bq (Σ) + b2s+q (Σ)).
2π
(B.1.20)
The remaining steps remains the same as the derivation for aq (Σ).

B.1.4 Proof for aq (Σi,j )

For the multivariate cases, we will index (B.1.2) with j, k and use the polari-
sation identity

1
hdp j , dpk it = hdp j + dpk it − hdp j it − hdpk it .


2
We first confirm that the polarization recovers the results we desire. We have

Z 2π
1
ak (dp j + dpk ) = cos(kt)(σ j (t) + σk (t))dW,
π 0

and

Z 2π
j k j 1 k
σ jj (t) + σkk (t) + 2σ jk (t) cos(kt) cos(lt)dt.


E[ ak (dp + dp ) al (dp + dp )] = 2
π 0

Therefore,

1h i
E[ ak (dp j + dpk ) al (dp j + dpk )] − E[ ak (dp j ) al (dp j )] − E[ ak (dpk ) al (dpk )]
2
Z 2π
1
= 2 σ jk (t) cos(kt) cos(lt)dt.
π 0
(B.1.21)
We note the linearity of Fourier transforms and thus ak (dp j + dpk ) = ak (dp j ) +
ak (dpk ), which means

ak (dp j + dpk ) al (dp j + dpk )

= ak (dp j ) + ak (dpk ) al (dp j ) + al (dpk )
  

= ak (dp j ) al (dp j ) + ak (dp j ) al (dpk ) + ak (dpk ) al (dp j ) + ak (dpk ) al (dpk ).

(B.1.22)
Using (B.1.22), (B.1.21) can be simplified into
62 Appendix B. Appendix Derivation

1h j k j k j j k k
i
E[ ak (dp + dp ) al (dp + dp )] − E[ ak (dp ) al (dp )] − E[ ak (dp ) al (dp )]
2
1h j k k j
i
= E[ ak (dp ) al (dp )] + E[ ak (dp ) al (dp )]
2
Z 2π
1
= 2 σ jk (t) cos(kt) cos(lt)dt
π 0
1
= ( a|k−l | (Σ jk ) + ak+l (Σ jk )).
2π
(B.1.23)
The 2π in the last equation of (B.1.23) is the 2π in the numerator of aq (Σi,j )
which was left out in [1].

B.1.5 Proof for a0 (Σi,j )

For a0 (Σi,j ), we look only at ( a2s (dp) + bs2 (dp)) and using the polarisation
identity

1h 2 j k 2 j k 2 j 2 j 2 k 2 k
i
E[ as (dp + dp ) + bs (dp + dp )] − E[ as (dp ) + bs (dp )] − E[ as (dp ) + bs (dp )] .
2
(B.1.24)
Using the linearity of Fourier transforms, the first expectation in (B.1.24) be-
comes

E as (dp j ) + as (dpk ) as (dp j ) + as (dpk ) + bs (dp j ) + bs (dpk ) bs (dp j ) + bs (dpk )







= E[ a2s (dp j )] + E[ a2s (dpk )] + 2E[ as (dp j ) as (dpk )]

+ E[bs2 (dp j )] + E[bs2 (dpk )] + 2E[bs (dp j )bs (dpk )].
(B.1.25)
Using this in (B.1.24), we get the polarisation of ( a2s (dp) + bs2 (dp)) as

1
2E[ as (dp j ) as (dpk ) + bs (dp j )bs (dpk )]. (B.1.26)
2

Therefore the scaling factors for a0 (Σ) and a0 (Σi,j ) are the same.

B.1.6 Proof for bq (Σi,j )

We have

Z 2π
1
ak (dp j + dpk ) = cos(kt)(σ j (t) + σk (t))dW,
π 0
B.2. Proof for Theorem 2.2.1 63

Z 2π
j k 1
bl (dp + dp ) = sin(lt)(σ j (t) + σk (t))dW,
π 0

and

Z 2π
j k j 1k
σ jj (t) + σkk (t) + 2σ jk (t) cos(kt) sin(lt)dt.


E[ ak (dp + dp )bl (dp + dp )] = 2
π 0

Therefore

1h j k j k j j k k
i
E[ ak (dp + dp )bl (dp + dp )] − E[ ak (dp )bl (dp )] − E[ ak (dp )bl (dp )]
2
Z 2π
1
= 2 σ jk (t) cos(kt) sin(lt)dt.
π 0
(B.1.27)
Using the linearity of Fourier transforms, we can simplify (B.1.27) to

1h i 1 2π
Z
E[ ak (dp j )bl (dpk )] + E[ ak (dpk )bl (dp j )] = 2 σ jk (t) cos(kt) sin(lt)dt
2 π 0
1
= (b (Σ jk ) + bk+l (Σ jk )).
2π |k−l |
(B.1.28)
We note that the scaling factors we recovered are the same scaling factors
recovered by [16].

B.2 Proof for Theorem 2.2.1

We present our own version for the proof that is different to [2], but the meth-
ods from [1] are used.
Assume first A-II and for ease of derivation that the price process follows a
one-dimensional Itô process

dp(t) = β(t)dt + σ (t)dW. (B.2.1)

Furthermore, we can assume that the drift term has no contribution which is
proven in [2]. Therefore the price process is given by

dp(t) = σ(t)dW. (B.2.2)

Further define
64 Appendix B. Appendix Derivation

Z 2π
1
ck (dp) = e−ikt dp(t).
2π 0

Now using the Itô isometry

1 h Z 2π Z 2π i
−ilt −ikt
E(cl ck ) = E e dp ( t ) e dp ( t )
(2π )2 0 0
1 h Z 2π Z 2π i
−ilt −ikt
= E e σ ( t ) e σ ( t )
(2π )2 0 0
(B.2.3)
Z 2π
1 h
−i ( l + k ) t
i
= E e Σ(t)dt
(2π )2 0
1
= c ( Σ ).
2π l +k
q
For q ∈ Z, define the Random Variable UN as the discrete convolution of the
Fourier coefficients for the price process.

N
1
∑ ck ck−q .
q
UN := (B.2.4)
2N + 1 s=− N

Therefore, it follows that

N
1
∑
 q  
E UN = E ck ck−q
2N + 1 s=− N
(B.2.5)
1 1
= ∑
2N + 1 |k|≤ N 2π
c q ( Σ ),

and thus

 q
2πE UN = cq (Σ). (B.2.6)

By setting the drift component β(t) ≡ 0, we have that ck (dp) is a complex

martingale. Therefore by introducing R2 valued gaussian variables

Gk := ck (dp),

where E[ Gk ] = 0. We can employ the same arguments used in the derivation

of aq (Σ) and further noting that energy identity is

kΣk2L2 = ∑ ck (Σ) .

2
k

From which we get that

B.3. Proof for Theorem 2.3.1 65

h
q
2
i 1 h i 1
∑0
q
E UN − E (U N ) = E(c2k ) E(c2q−k0 ) + 2E(ck cq−k0 )2 − c q ( Σ )2
(2N + 1)2 k,k
(2π ) 2

q
2 1 1 1
= E UN +
(2N + 1)2 ∑0 2π ck+q+k0 (Σ) − (2π )2 cq (Σ)2
k,k
1
≤ kΣk2L2 .
2N + 1
(B.2.7)
By combining (B.2.6) and (B.2.7), we complete the theorem as

N
1 1
c (Σ) = lim
2π k ∑ cs (dp)ck−s (dp),
N →∞ 2N + 1 s=− N
(B.2.8)

where the convergence of the Bohr convolution is attained in probability.

B.3 Proof for Theorem 2.3.1

The proof is based on [3], we add a few explanatory steps and details to their
derivation
(i) We first assume that µl ≡ 0 and 0 ≤ t ≤ T. We now show that Un → θ in
L2 as n → ∞. We introduce some auxiliary symbols: let Kij := I{ I i ∩ J j 6=∅} . Fur-
thermore, for each measurable set I on [0, ∞), we define (signed) measures
by
Z
0
υ( I ) := υ ( I ) := σ1 σ2 ρdt,
Z I
k
υ ( I ) := (σk )2 dt, k = 1, 2.
I

We further introduce identities for k = 0, 1, 2

   
∑ I i 1{ I i 6 = ∅ } =
v k
∑ I i = vk ((0, T ]),
v k
i i
 
∑ vk I i ∩ J j Kij = vk ((0, T ]),
i,j
   
∑v k
I ∩Ji j
Kij = v I i , k
j

these identities hold due to the fact that I i and J j are the partitions over (0, T ].
Finally, for each measurable set I on [0, ∞), define

Z T
k
∆P ( I ) := 1 I (t)σk dW k , k = 1, 2.
0
66 Appendix B. Appendix Derivation

We can show that Un is unbiased, for each n,

" # " #
    o n  
E [Un ] = E ∑ E ∆Pl I i ∆P2 J j |Π Kij = E ∑ v I i ∩ J j Kij = θ.
i,j i,j

Remark B.3.1 We note that the inner expectation is for notation indicating that we
know Π and that the outer expectation can be brought in due to linearity. This is
to allow later parts of the derivation become more simple. The reason why we say
the inner expectation is for notation is because we can recover θ with just the outer
expectation:
" #
   
E[Un ] = E ∑ ∆Pl I i ∆P2 J j Kij
i,j
T Z T
Z 
= ∑E 1 I i (t)σ dW 1 1 2 2
1 J j (t)σ dW Kij
i,j 0 0
T
Z 
= ∑E 12
1 I i ∩ J j (t)σ dt Kij (B.3.1)
i,j 0
Z 
= ∑E 1 2
σ σ ρdt Kij
i,j Ii ∩ J j
" #
 
=E ∑v I i ∩ J j Kij = θ.
i,j

The third equation follows from the second in (B.3.1) using Itô’s Isometry. For the
remainder of the proof we adopt inner expectation used by [3] to make things simpler.
We now want to show that E Un2 = θ 2 + o (1). This would mean that
 

var [Un ] = o (1) so that Un → θ in L2 as n → ∞. To do so, we start with

" #
h i n     0  0 o
E Un2 =E ∑ i j i j
E ∆P I ∆P J ∆P I ∆P J |Π Kij Ki0 j0 ,
2 ll 2

i,j,i0 ,j0

for which the summations can be decomposed into

∑ = ∑ + ∑ + ∑ + ∑ =: D1 + D2 + D3 + D4 .
i,j,i0 ,j0 0 0 0 0 0 0 0 0
i, j, i , j : i, j, i , j : i, j, i , j : i, j, i , j :
i0 = i, j0 = j i0 = i, j0 6= j i0 6= i, j0 = j i0 6= i, j0 6= j

We will use this decomposition to calculate the expectation using four cases.
B.3. Proof for Theorem 2.3.1 67

Case 1: i = i0 , j = j0 . Let L1 := I i ∩ J j , L2 := I i \ L1 and L3 := J j \ L1 . We have


 2  2 
D1 = ∑ E ∆P I i
∆P J j |Π Kij .
2 1
i,j

Using the independence of the increments and identities such as vk ( I i ) =

vk ( L2 ) + vk ( L1 ) and vk ( J j ) = vk ( L3 ) + vk ( L1 ), we get that

  2  2   2  2 
i j
E ∆P I l
∆P J
2
|Π = E ∆P ( L2 ) + ∆P ( L1 )
1 1
∆P ( L3 ) + ∆P ( L1 ) |Π
2 2

n o n o
2 2 2 2
= E ∆P ( L2 ) ∆P ( L1 ) |Π + E ∆P ( L1 ) ∆P ( L1 ) |Π
1 2 1 2
n o n o
+E ∆P1 ( L1 )2 ∆P2 ( L3 )2 |Π + E ∆P1 ( L2 )2 ∆P2 ( L3 )2 |Π
= v1 ( L2 ) v2 ( L1 ) + 2v ( L1 )2 + v1 ( L1 ) v2 ( L1 )
+ v1 ( L1 ) v2 ( L3 ) + v1 ( L2 ) v2 ( L3 )
= [v1 ( L2 ) + v1 ( L1 )]v2 ( L1 ) + [v1 ( L1 ) + v1 ( L2 )]v2 ( L3 ) + 2v( L1 )2
= [v1 ( L1 ) + v1 ( L2 )][v2 ( L1 ) + v2 ( L3 )] + 2v( L1 )2
= v1 ( I i )v2 ( J j ) + 2v( I i ∩ J j )2 .

The first equation comes from the definition that I i = L2 + L1 and J j = L3 +

L1 . The second equation comes from the fact that the cross products from
each of the quadratic terms are 0. The third equation follows from the fact
that for any (deterministic) interval I, ∆P1 ( I ) and ∆P2 ( I ) are jointly normal
with respective mean and variance 0 and vk ( I ), k = 1, 2, and with covariance
v( I ). Thus E[∆P1 ( I )∆P2 ( I )] = 2v( I )2 + v1 ( I )v2 ( I ) using some Multivariate
reults. Therefore we have,

     2
D1 = ∑ v I v J Kij + 2 ∑ v I ∩ J
i12 j i j
Kij . (B.3.2)
i,j i,j

Now looking at the first term on the right hand side of (B.3.2) and noting that
the σk are bounded, we get

    Z  2  Z  2 
∑ I i v2 J j Kij =
v 1
∑ Ii
σ 1
dt
Jj
σ 2
dt Kij
i,j i,j
 2  2
≤ sup σ 1
sup σ 2
∑ | I i || J j |Kij .
0≤ t ≤ T 0≤ t ≤ T i,j

We now want to show that

68 Appendix B. Appendix Derivation

E ∑ | I i || J j |Kij = o (1).
i,j

To do so, we decompose

∑ | I i || J j |Kij = ∑ | I i || J j |Kij 1{| I i |≥| J j |} + ∑ | I i || J j |Kij 1{| I i |<| J j |} .

i,j i,j i,j

Noting that ∑ j | J j |Kij 1{| I i |≥| J j |} ≤ 3| I i | for some fixed i, we have

∑ | I i || J j |Kij 1{| I i |≥| J j |} = ∑ | I i | ∑ | J j |Kij 1{| I i |≥| J j |} ≤ 3 ∑ | I i |2,

i,j i j i

hence,

E ∑ | I i || J j |Kij 1{| I |≥| J |} ≤ 3E ∑ | I i |2 .

i,j i

By symmetry, we have

E ∑ | I i || J j |Kij ≤ 3E ∑ | I i |2 + 3E ∑ | J j |2 . (B.3.3)
i,j i j

We see that (B.3.3) is o (1) under Condition (A-IV) (ii), (C(ii) in remark 3.1 of
[3]). Similarly, we can ascertain that for any random partition ( Ĩ i ) of (0, T ]
satisfying (A-IV) (ii),

 2
E ∑ v Ĩ i = o (1). (B.3.4)
i

Thus the second term on the right hand side of (B.3.2) can be shown to be
of o P (1) by choosing ( I i ∩ J j ) as the partition. Hence it follows that E[ D1 ] =
o (1).
Case 2: i = i0 , j 6= j0 . This yields
  2    0 
D2 = ∑ E ∆P I il
∆P J ∆P J j |Π Kij Kij0 .
2 j 2

i,j0 :j6= j0

0
Let L1 := I i ∩ J j , L2 := I i ∩ J j , and L3 := I i \ ( L1 ∪ L2 ). Then using the
independence of increments,
B.3. Proof for Theorem 2.3.1 69

  2    0 
i
E ∆P I l
∆P J ∆P J j |Π
2 j 2

  2 
i
= E ∆P I ∆P ( L1 ) ∆P ( L2 ) |Π
l 2 2

 2 
=E ∆P ( L1 ) + ∆P ( L3 ) + ∆P ( L2 ) ∆P ( L1 ) ∆P ( L2 ) |Π
1 1 1 2 2

n o n o
= 2E ∆P1 ( L1 ) ∆P2 ( L1 ) |Π E ∆Pl ( L2 ) ∆P2 ( L2 ) |Π
   
i j i j0
= 2v ( L1 ) v ( L2 ) = 2v I ∩ J v I ∩ J .

The third equation follows because when we expand the second equation,
all the non-overlapping increments reduces to 0 using the independence of
increments and that for any (deterministic) interval I, ∆P1 ( I ) and ∆P2 ( I ) are
jointly normal with respective mean and variance 0 and vk ( I ), k = 1, 2, and
with covariance v( I ).
Hence,
   
j0
D2 = 2 ∑ i i j
v I ∩ J v I ∩ J Kij Kij0
i,j,j0 :j0 6= j
( !)
     
j0
= 2∑ ∑v I i ∩ J j Kij ∑0 v Ii ∩ J Kij0 − v I i ∩ J j
i j j
 2  2
= 2 ∑ v I i − 2 ∑ ∑ v I i ∩ J j Kij ,
i i j

which follows becuase of the relation ∑ j vk I i ∩ J j Kij = vk I i . Finally, we

see that E[ D2 ] = o (1) by using (B.3.4) and the fact that ( I i ∩ J j ) partitions
]0, T ].
Case 3: i 6= i0 , j = j0 . The same argument in case 2 applies here by symmetry,
thus we can obtain E[ D3 ] = o (1).
0 0
Case 4: i 6= i0 , j 6= j0 . Let L1 := I i ∩ J j , L2 := I i ∩ J j . Note that for i, i0 , j, j0
such that i 6= i0 , j 6= j0 and Kij Ki0 j0 = 1 means that Ki0 j Kij0 = 0. Furthermore,
due to the identity
  
1 − Ki 0 j 1 − Kij0 + Ki0 j + Kij0 ≡ 1,
n 0 0
o
we can decompose the event {Kij Ki0 j0 = 1} further into three subcases, ∩ = ∅, Ii Jj Ii ∩ Jj =∅ ,
n 0 0
o n 0 0
o
I i ∩ J j 6= ∅, I i ∩ J j = ∅ and I i ∩ J j = ∅, I i ∩ J j 6= ∅ each of which re-
n   o n o n o
spectively corresponds to 1 − Ki0 j 1 − Kij0 = 1 , Ki0 j = 1 and Kij0 = 1 .
n 0 0
o
Case 4(a): I i ∩ J j = ∅, I i ∩ J j = ∅ . We have by analogy with case 2
70 Appendix B. Appendix Derivation

n      0  0 o   
∑ E ∆P1 I i ∆P2 J j ∆P1 I i ∆P2 J j |Π Kij Ki0 j0 1 − Ki0 j 1 − Kij0
i,j,i0 ,j0 :i 6=i0 ,j6= j0
n o   
= ∑ E ∆P ( L1 ) ∆P ( L1 ) ∆P ( L2 ) ∆P ( L2 ) |Π Kij Ki0 j0 1 − Ki0 j
1 2 1 2
1 − Kij0
i,j,i0 ,j0 :i 6=i0 ,j6= j0
  
= ∑ v( L1 )v( L2 )Kij Ki0 j0 1 − Ki0 j 1 − Kij0
i,j,i0 ,j0 :i 6=i0 ,j6= j0

n 0 0
o 0
Case 4(b): I i ∩ J j 6= ∅, I i ∩ J j = ∅ . Let L3 := I i ∩ J j , L4 := J j \ ( L1 ∪ L3 )
0
and L5 := I i \ ( L2 ∪ L3 ),

E ∆P1( L1 ) ∆P2 ( L1 ) + ∆P2 ( L4 ) + ∆P2 ( L3 ) ∆P1 ( L3 ) + ∆P1 ( L5 ) + ∆P1 ( L2 ) ∆P2 ( L2 ) |Π



= E ∆P1 ( L1 ) ∆P2 ( L1 ) ∆P1 ( L2 ) ∆P2 ( L2 ) |Π = v ( L1 ) v ( L2 )

Therefore,

n      0  0 o
∑ E ∆P1 I i ∆P2 J j ∆P1 I i ∆P2 J j |Π Kij Ki0 j0 Ki0 j
i,j,i0 ,j0 :i 6=i0 ,j6= j0

= ∑ v ( L1 ) v ( L2 ) Kij Ki0 j0 Ki0 j .

i,j,i0 ,j0 :i 6=i0 ,j6= j0

n 0 0
o
Case 4(c): ∩ = ∅, Ii Jj Ii ∩ Jj 6= ∅ . By symmetry, we can obtain using the
same technique as 4(b):

∑ v ( L1 ) v ( L2 ) Kij Ki0 j0 Kij0 .

i,j,i0 ,j0 :i 6=i0 ,j6= j0

Combining all three subcases together, we get

h ih    i
D4 = ∑ v ( L1 ) v ( L2 ) Kij Ki0 j0 1 − Ki 0 j 1 − Kij0 + Kij0 + Kij0
i,j,i0 ,j0 :i 6=i0 ,j6= j0
  0 0
= ∑0 v I i ∩ J j v( I i ∩ J j )Kij Ki0 j0
i,j,i0 ,j :i 6=i ,j6= j0
0
   
i0 j0
= ∑v I ∩ J i j
Kij ∑ v ( I ∩ J ) Ki 0 j 0 .
i,j i0 ,j0 :i0 6=i,j0 6= j

Now for fixed i and j,

B.3. Proof for Theorem 2.3.1 71

0 0
∑ v ( I i ∩ J j ) Ki 0 j 0
i0 ,j0 :i 6=i,j0 6= j
0
 0 0
    0
  0 
= ∑ v I i ∩ J j Ki 0 j 0 − v I i ∩ J j − ∑ v I i ∩ J j Kij0 − ∑ v I i ∩ J j Ki 0 j
i0 ,j0 j0 :j0 6= j ii :i0 6=i
 0 0
  
= ∑ v I i ∩ J j Ki 0 j 0 − v I i ∩ J j − v ( I i ) − v ( J j ),
i0 ,j0

we note the typo in [3], they have v I i ∩ J j instead of −v I i ∩ J j . We now

have

   0   2
j0
D4 = ∑ v I i
∩ J j
K ij ∑ v I i
∩ J K 0
ij 0 − ∑ v I i
∩ J j
Kij
i,j i0 ,j0 i,j
   
− ∑ v I i ∩ J j Kij v( I i ) − ∑ v I i ∩ J j Kij v( J j )
i,j i,j
 2
= v(]0, T ])2 − ∑ v I i ∩ J j Kij − ∑ v( I i )2 − ∑ v( J j )2 .
i,j i j

Thus E[ D4 ] = θ 2 + o (1) by (B.3.4).

Therefore, E[Un2 ] = E[ D1 + D2 + D3 + D4 ] = θ 2 + o (1).
(ii) Now we consider the case with non-zero drift such that sup0≤t≤T |uk | ∈
R. R.
L4 , k = 1, 2. Let A.k := 0 uk dt, M.k := 0 σk dW k , k = 1, 2 and

       
B0 := ∑ ∆M I ∆M J j Kij ,
i 2 1
B1 := ∑ ∆A I ∆M J j Kij ;
i 2 1
i,j i,j
       
B2 := ∑ ∆M1 I i ∆A2 J j Kij , B3 := ∑ ∆A 1
I i
∆A 2
J j Kij .
i,j i,j

Note that

Z Z
! Z Z
| B1 | = ∑ Ii
µ1 dt ∑ Jj
σ2 dW 2 Kij ≤∑
Ii
µ1 dt ∑ Jj
σ2 dW 2 Kij
i j i j
t
Z 
1 2 2 i j
≤ T sup µ · max sup σ dW |t − s| ≤ I + 2 max J , s, t ∈ [0, T ]
0≤ t ≤ T i s j
Z t 
1 2 2 i j
≤ T sup µ · sup σ dW |t − s| ≤ max I + 2 max J , s, t ∈ [0, T ] ,
0≤ t ≤ T s i j
(B.3.5)
72 Appendix B. Appendix Derivation

from which we see that B1 is in L2 because σ2 is bounded and the supremum

of µ1 is in L4 . Moreover, under Condition (A-IV) (ii), E[ B12 ] = o (1) as n → ∞
by the dominated convergence theorem. E[ B22 ] = o (1) and E[ B32 ] = o (1) can
be shown similarly.
Furthermore,
h i h i h i
E (Un − θ )2 ≤ 2E ( B0 − θ )2 + 8E B12 + B22 + B32 ,

and together with (i), we have shown that Un → θ in L2 as n → ∞ because

B0 is just the driftless case.
Showing Un → θ in probability as n → ∞ follows from (B.3.5), as all that is
needed is Bk → 0 in probability as n → ∞, k=1,2,3.
 73

Appendix C

Miscellaneous

C.1 Run times and General issues

We compare first the various run times of algorithm 1, 2 and 3. We begin
the experiment by first simulating 500 synchronous seconds from a bivari-
ate GBM using algorithm 10, we then sample individual samples paths with
an exponential inter-arrival time with a mean of 15 seconds to obtain asyn-
chronous sample paths with an average of 33 data points. With these sample
paths we measure the run-times for the various algorithms 1 100 times and
the distribution for each algorithm is shown in figure C.1.

(a) Synchronous Run Time (b) Asynchronous Run Time

HY HY

MM Complex MM Complex

MM Trig MM Trig

100 1000 10000 3 10 30 100

Time [milliseconds] Time [milliseconds]

F IGURE C.1: We compare the run times of algorithm 1, 2 and 3.

Specifically, (a) looks at the run time for 500 synchronous data
points while (b) looks at the run time for an average of 33 data
points where the asynchrony is induced by exponential inter-
arrival time samples of the synchronous data points in (a). It
is clear that the Complex MM has the fastest run time while
HY has a significantly slower run time compared to the MM
algorithms.
1 Forthe computation of the MM estimators, the number of Fourier coefficients is chosen
to be the Nyquist frequency based on the highest available sampling frequency present in
the sample.
74 Appendix C. Miscellaneous

From figure C.1 2 , we see that for both the synchronous case and the asyn-
chronous case, the HY algorithm is order of magnitudes longer than the MM
algorithms. This is due to the Kanatani Weight matrix using a double for-
loop to check for overlapping intervals. Furthermore, we notice that algo-
rithm 2 is faster than algorithm 1, this is because algorithm 2 has been sup-
plemented with Rcpp and RcppArmadillo code to improve the computation
time.
We started experimenting with C++ code because memory issues started oc-
curring with base R when computing the empirical data. This was because all
the Fourier coefficients are computed with one matrix multiplication for com-
putation efficiency; however, due to the large nature of empirical data, along
with the fact that we computed Fourier coefficients based on the Nyquist
frequency for the highest available sampling frequency present in the data -
computing one pair of stock for one day of data required an matrix of dimen-
sions [15,000 x 1,500,000] to be initialised.
Upon further investigation, we found that R uses 8 bytes of memory to store a
double precision float, therefore a matrix with dimensions [15,000 x 1,500,000]
uses 167.6Gb of memory to store the object - meaning that packages such as
bigmemory which allows the data structure to be allocated to shared mem-
ory was not going to be particularly helpful due to physical constraints of
our hardware. Thus the only solution left was to remove the vectorisation
and compute each Fourier coefficients using a single for-loop.
It must be noted, even though the vectorisation was removed from algorithm
2, it took 4 days to compute the correlation matrix for 1 week of data of the
10 assets considered, while algorithm 3 took 7 days. Furthermore, it must be
noted that algorithm 3 will have memory issues that are not easily solvable as
the Kanatani weight matrix becomes larger when considering longer periods
such as the correlation for a month of data.
These problems are non-trivial computer science problems which severely
impact the real-time implementation of these two estimators. Fixing these
problems either requires very efficient parallelisation [5] or more efficient al-
gorithms are needed. More effective algorithms for 2 include a Fast Fourier
Transform (FFT) for the synchronous case and a Non-Uniform FFT (NUFFT)
for the more practical asynchronous case.
Remark C.1.1 One important point to notice is that we did not use algorithm 1
to compute any correlation results. This is because the algorithm is producing the
wrong correlation estimates for the asynchronous cases. The exact reason for this
error has not yet been figured out. We know that it is algorithm 1 which is wrong
rather than algorithm 2 because for the asynchronous case, the integrated variance
should be the same as the HY estimate (since the variance is synchronous) - however
it is not.
2 Figure C.1 can be reproduced using Compute Time.R
C.2. Supporting plots 75

C.2 Supporting plots

Geometric Brownian Motion Merton Model
103
102.0 0.004

101.5 102
Price (P)

Price (P)
101.0 0.0025
101
100.5 0.002

Sample (Merton)
Sample (GBM)
100.0 100

0 100 200 300 400 500 0 100 200 300 400 500
Time Time 0.0000

0.000
0.004

0.0025
Return(log∆P)

Return(log∆P)
0.002

0.0000 −0.0025
0.000 −0.002

−0.002 −0.0025

0 100 200 300 400 500 −2 0 2 0 100 200 300 400 500 −2 0 2
Time Theoretical N(0,1) Time Theoretical N(0,1)

( A ) Geometric Brownian motion ( B ) Merton Model

GARCH (1,1) − Reno GARCH (1,1) − Andersen

100 100
Price (P)

Price (P)
98 98
0.005 0.005

96 96
Sample (Garch(1,1))

Sample (Garch(1,1))
0 100 200 300 400 500 0 100 200 300 400 500
Time 0.000 Time 0.000

0.005 0.005
Return(log∆P)

Return(log∆P)

−0.005 −0.005
0.000 0.000

−0.005 −0.005

0 100 200 300 400 500 −2 0 2 0 100 200 300 400 500 −2 0 2
Time Theoretical N(0,1) Time Theoretical N(0,1)

( C ) GARCH(1.1) specification of [12] ( D ) GARCH(1.1) specification of [23]

Ornstein Uhlenbeck Variance Gamma
100.015
4e−05
100.5 0.004
100.010
Price (P)

Price (P)

100.0 100.005

100.000
Sample (Ornstein Uhlenbeck)

2e−05
Sample (Variance Gamma)

99.5 0.002
99.995

99.0 99.990
0 100 200 300 400 500 0 100 200 300 400 500
Time Time
0.000 0e+00
0.004 4e−05
Return(log∆P)

Return(log∆P)

0.002 2e−05

0.000 −0.002 0e+00 −2e−05

−0.002 −2e−05

0 100 200 300 400 500 −2 0 2 0 100 200 300 400 500 −2 0 2
Time Theoretical N(0,1) Time Theoretical N(0,1)

( E ) Ornstein Uhlenbeck ( F ) Variance Gamma

F IGURE C.2: Price paths, returns and QQ-plot for the various
SDEs.

Figure C.2 highlights the resulting price paths generated using the Algo-
rithms from appendix A, and further shows the returns and QQ-plots as-
sociated with the various SDEs considered in this paper.
76 Appendix C. Miscellaneous

(a) MM 1 HR Close Price (e) HY 1 HR Close Price

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.1411 −0.5 AGL |ρij| = 0.1413 −0.5

NPN NPN
BTI BTI
−1.0 −1.0

MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI
(b) MM 10 Min Close Price (f) HY 10 Min Close Price
ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.0734 −0.5 AGL |ρij| = 0.08 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

(c) MM 1 Min Close Price (g) HY 1 Min Close Price

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.0379 −0.5 AGL |ρij| = 0.0469 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

(d) MM TAQ (h) HY TAQ

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.0355 −0.5 AGL |ρij| = 0.1051 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

F IGURE C.3: Demonstrating the existence of the Epps effect

when returns are all computed on the TAQ data rather than on
the closing bar prices. From (a) to (d), we have the MM estima-
tor applied to TAQ returns sampled every 1 hour, 10 minute,
1 minute and all TAQ return samples respectively using algo-
rithm 2. From (e) to (h), we have the HY estimator applied to
TAQ returns sampled every 1 hour, 10 minute, 1 minute and
all the TAQ return samples respectively using algorithm 3. The
Epps effect is still present even when all the returns are com-
puted based on the highest available sampling frequency rather
than on the respective sampling intervals.
C.2. Supporting plots 77

(a) MM 1 HR VWAP (e) HY 1 HR VWAP

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.1418 −0.5 AGL |ρij| = 0.1449 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI
(b) MM 10 Min VWAP (f) HY 10 Min VWAP
ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.097 −0.5 AGL |ρij| = 0.0961 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

(c) MM 1 Min VWAP (g) HY 1 Min VWAP

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.0627 −0.5 AGL |ρij| = 0.0784 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

(d) MM TAQ (h) HY TAQ

ρ ρ
FSR FSR
1.0 1.0

SHP SHP
ABG ABG
0.5 0.5
NED NED
SBK SBK
0.0 0.0
SOL SOL
MNP MNP
AGL |ρij| = 0.0355 −0.5 AGL |ρij| = 0.1051 −0.5

NPN NPN
BTI BTI
−1.0 −1.0
MNP

ABG

MNP

ABG
NPN

SBK
NED

NPN

SBK
NED
AGL

SOL

SHP

AGL

SOL

SHP
FSR

FSR
BTI

BTI

F IGURE C.4: Demonstrating the existence of the Epps effect

when returns are all computed on the TAQ data and the VWAP
aggregation is performed on the TAQ returns. From (a) to (d),
we have the MM estimator applied to the VWAP aggregation
of TAQ returns from 1 hour, 10 minute, 1 minute bars and pure
TAQ samples respectively using algorithm 2. From (e) to (h),
we have the HY estimator applied to the VWAP aggregation of
TAQ returns from 1 hour, 10 minute, 1 minute bars and pure
TAQ samples respectively using algorithm 3. The Epps effect
is still present even under this bizarre aggregation of TAQ re-
turns.
78 Appendix C. Miscellaneous

Figures C.3 and C.4 show the initial mistake we made with figures 5.1 and
5.2 respectively. The mistake was computing the returns before computing
the OHCLV prices using algorithm 5, therefore we initially inputted the re-
turns computed on the TAQ data into algorithm 5. Therefore the bar data in
figure C.3 is rather the TAQ return samples sampled every 1 hour, 10 minute
and 1 minute respectively while the bar data in figure C.4 is the VWAP algo-
rithm applied to the TAQ return samples for 1 hour, 10 minute and 1 minute
intervals respectively.
What is interesting about this mistake is that it still shows the existence of
the Epps effect even though all the returns are computed based on the high-
est available sampling frequency rather than ranging sampling intervals. The
only commonality with the samples in figures C.3, C.4 and 5.1, 5.2 is that the
number of samples decrease as the sampling intervals decrease. This begs the
question regarding the relationship between the Epps effect and the number
of samples considered. Additionally, are there any factors other than asyn-
chrony, lead-lag and smaller sampling intervals which contribute towards
the Epps effect? Unfortunately, we have not been able to answer these ques-
tions.
 79

Bibliography

[1] P. Malliavin and M. E. Mancino, “Fourier series method for measure-

ment of multivariate volatilities”, Finance and Stochastics, vol. 6, no. 1,
pp. 49–61, 2002, ISSN: 0949-2984. DOI: 10 . 1007 / s780 - 002 - 8400 - 6.
[Online]. Available: [Link]
[2] ——, “A fourier transform method for nonparametric estimation of
multivariate volatility”, Ann. Statist., vol. 37, no. 4, pp. 1983–2010, Aug.
2009. DOI: 10.1214/08-AOS633. [Online]. Available: [Link]
10.1214/08-AOS633.
[3] T. Hayashi and N. Yoshida, “On covariance estimation of non-synchronously
observed diffusion processes”, Bernoulli, vol. 11, no. 2, pp. 359–379,
Apr. 2005. DOI: 10.3150/bj/1116340299. [Online]. Available: https:
//[Link]/10.3150/bj/1116340299.
[4] T. W. Epps, “Comovements in stock prices in the very short run”, Jour-
nal of the American Statistical Association, vol. 74, no. 366, pp. 291–298,
1979, ISSN: 01621459. [Online]. Available: http : / / www . jstor . org /
stable/2286325.
[5] D. Hendricks, “Using real-time cluster configurations of streaming asyn-
chronous features as online state descriptors in financial markets”, Pat-
tern Recognition Letters, vol. 97, pp. 21 –28, 2017, ISSN: 0167-8655. DOI:
https : / / doi . org / 10 . 1016 / j . patrec . 2017 . 06 . 026. [Online].
Available: [Link]
S0167865517302325.
[6] I. Mastromatteo, M. Marsili, and P. Zoi, “Financial correlations at ultra-
high frequency: Theoretical models and empirical estimation”, The Eu-
ropean Physical Journal B, vol. 80, no. 2, pp. 243–253, 2011, ISSN: 1434-
6036. DOI: 10.1140/epjb/e2011-10865-y. [Online]. Available: https:
//[Link]/10.1140/epjb/e2011-10865-y.
[7] E. Barucci and R. Renò, “On measuring volatility and the garch fore-
casting performance”, Journal of International Financial Markets, Institu-
tions and Money, vol. 12, pp. 183–200, Jul. 2002. DOI: 10.1016/S1042-
4431(02)00002-1.
[8] M. Mancino, M. Recchioni, and S. Sanfelici, Fourier-Malliavin Volatility
Estimation Theory and Practice. Springer International Publishing, Mar.
2017, ISBN: 978-3-319-50967-9. DOI: 10.1007/978-3-319-50969-3.
[9] M. Mancino and S. Sanfelici, “Robustness of fourier estimator of in-
tegrated volatility in the presence of microstructure noise”, Computa-
tional Statistics Data Analysis, vol. 52, pp. 2966–2989, Feb. 2008. DOI:
10.1016/[Link].2007.07.014.
80 Bibliography

[10] J. E. Griffin and R. C. Oomen, “Covariance measurement in the pres-

ence of non-synchronous trading and market microstructure noise”,
Journal of Econometrics, vol. 160, no. 1, pp. 58–68, 2011. [Online]. Avail-
able: https : / / ideas . repec . org / a / eee / econom / v160y2011i1p58 -
[Link].
[11] Y. AÃt-Sahalia, J. Fan, and D. Xiu, “High-frequency covariance esti-
mates with noisy and asynchronous financial data”, Journal of the Amer-
ican Statistical Association, vol. 105, pp. 1504–1517, Jun. 2010. DOI: 10.
2139/ssrn.1631344.
[12] R. Renò, “A closer look at the epps effect”, International Journal of The-
oretical and Applied Finance, vol. 06, Nov. 2001. DOI: 10 . 2139 / ssrn .
314723.
[13] O. V. Precup and G. Iori, “Cross-correlation measures in the high-frequency
domain”, The European Journal of Finance, vol. 13, no. 4, pp. 319–331,
2007. DOI: 10.1080/13518470600813565. eprint: [Link]
1080/13518470600813565. [Online]. Available: [Link]
1080/13518470600813565.
[14] E. Derman, “The perception of time, risk and return during periods of
speculation”, Quantitative Finance, vol. 2, no. 4, pp. 282–296, 2002. DOI:
10.1088/1469- 7688/2/4/304. eprint: [Link]
1469-7688/2/4/304. [Online]. Available: [Link]
1469-7688/2/4/304.
[15] C. Malherbe, “Fourier method for the measurement of univariate and
multivariate volatility in the presence of high frequency data”, Mas-
ter’s thesis, University of Cape Town, 2007.
[16] T. Hoshikawa, K. Nagai, T. Kanatani, and Y. Nishiyama, “Nonparamet-
ric estimation methods of integrated multivariate volatilities”, Econo-
metric Reviews, vol. 27, no. 1-3, pp. 112–138, 2008. DOI: 10.1080/07474930701853855.
eprint: https : / / doi . org / 10 . 1080 / 07474930701853855. [Online].
Available: [Link]
[17] D. Hendricks, T. Gebbie, and D. Wilcox, “High-speed fourier method
estimation of covariances from asynchronous data”, working paper,
2017.
[18] T. Gebbie, D. Wilcox, C. Malherbe, and D. Hendricks, Ftcorrgpu.m, 2005.
[19] T. Kanatani, “Optimally weighted realized volatility”, Jan. 2004.
[20] B. Tóth and J. Kertész, “The epps effect revisited”, Quantitative Finance,
vol. 9, no. 7, pp. 793–802, 2009. DOI: 10 . 1080 / 14697680802595668.
eprint: https : / / doi . org / 10 . 1080 / 14697680802595668. [Online].
Available: [Link]
[21] W Moon and J. S. Wettlaufer, “On the interpretation of stratonovich
calculus”, New Journal of Physics, vol. 16, no. 5, p. 055 017, 2014. DOI:
10.1088/1367-2630/16/5/055017. [Online]. Available: [Link]
org/10.1088%2F1367-2630%2F16%2F5%2F055017.
[22] P. Glasserman, Monte Carlo methods in financial engineering. New York:
Springer, 2004, ISBN: 0387004513 9780387004518 1441918221 9781441918222.
[Online]. Available: [Link]
Bibliography 81

Stochastic-Modelling-Probability/dp/0387004513/ref=pd_sim_b_
68?ie=UTF8&refRID=1AN8JXSDGMEV2RPHFC2A.
[23] T. G. Andersen and T. Bollerslev, “Answering the skeptics: Yes, stan-
dard volatility models do provide accurate forecasts”, International Eco-
nomic Review, vol. 39, no. 4, pp. 885–905, 1998, ISSN: 00206598, 14682354.
[Online]. Available: [Link]
[24] T. G. Andersen and T. Teräsvirta, “Realized volatility”, in Handbook
of Financial Time Series, T. Mikosch, J.-P. Kreiß, R. A. Davis, and T. G.
Andersen, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2009,
pp. 555–575, ISBN: 978-3-540-71297-8. DOI: 10.1007/978-3-540-71297-
8_24. [Online]. Available: [Link] 3- 540-
71297-8_24.
[25] A. G. Hawkes, “Spectra of some self-exciting and mutually exciting
point processes”, Biometrika, vol. 58, no. 1, pp. 83–90, 1971, ISSN: 00063444.
[Online]. Available: [Link]
[26] R. Arbi, “A reproducible approach to equity backtesting”, Master’s the-
sis, University of Cape Town, 2019.
[27] Volume 00e-trading and information overview, Johannesburg Stock Ex-
change, Sandown, Johannesburg, South Africa, 2019.
[28] D. Easley, R. F. Engle, M. O’Hara, and L. Wu, “Time-Varying Arrival
Rates of Informed and Uninformed Trades”, Journal of Financial Econo-
metrics, vol. 6, no. 2, pp. 171–207, Feb. 2008, ISSN: 1479-8409. DOI: 10.
1093 / jjfinec / nbn003. eprint: http : / / oup . prod . sis . lan / jfec /
article - pdf / 6 / 2 / 171 / 2594496 / nbn003 . pdf. [Online]. Available:
[Link]
[29] D. Easley, M. M. López de Prado, and M. O’Hara, “The volume clock:
Insights into the high-frequency paradigm”, The Journal of Portfolio Man-
agement, vol. 39, no. 1, pp. 19–29, 2012, ISSN: 0095-4918. DOI: 10.3905/
jpm.2012.39.1.019. eprint: [Link]
39/1/[Link]. [Online]. Available: [Link]
com/content/39/1/19.
[30] D. Easley, M. Lopez de Prado, and M. O’Hara, “Flow toxicity and liq-
uidity in a high frequency world”, Review of Financial Studies, vol. 25,
Feb. 2012. DOI: 10.2139/ssrn.1695596.
[31] T. Gebbie, D. Wilcox, C. Malherbe, and D. Hendricks, Fftcorrgpu.m,
2005.
[32] P. Chang, R. Bukuru, and T. Gebbie, 2019. [Online]. Available: https:
//[Link]/rogerbukuru/Honours-Project.
[33] B. Oksendal, Stochastic Differential Equations (3rd Ed.): An Introduction
with Applications. Berlin, Heidelberg: Springer-Verlag, 1992, ISBN: 3-387-
53335-4.