Some Time Series LispTick examples
Returns
Compare two timeserie returns during the day.
Create a function ts asking for trade price for a stock for a day.
Create another function tsreturn computing return by dividing by open value and substarcting 1.
(defn ts[code start stop]
(timeserie @trade-price "refinitiv" code start stop))
(defn tsreturn[code start stop]
(- (/ (ts code start stop) (first (ts code start stop))) 1))
[
(graphsample 1920 (tsreturn "BNPP.PA" 2020-03-02T09:10 2020-03-02T17:30))
(graphsample 1920 (tsreturn "SOGN.PA" 2020-03-02T09:10 2020-03-02T17:30))
]
Interpolations
Use known values to interpolate for missing times.
By default LispTick considers timeserie as piecewise constant.
Here we present how to create intermediate points by using linear or cubic interpolations.
Those implementations ensure interpolation occurs only at requested points, i.e. for timestamp in timeserie named times.
Linear & Cubic
(def
start 2000-01-01
pi 3.1415926535
window 16s
steps 32
intersteps (* 32 steps))
;number of window from start
(defn sfs[t]
(/ (- t start) window))
;reference time serie cos(14πt²)+t²
;t = 0 at start and 1 at start + window
(def ts
(timeserie
(range start (+ start window) (/ window steps))
(fn[t] (+ (cos (* 7 2 pi (sqr (sfs t)))) (sqr (sfs t))))))
;desired timestamp for interpolation
(def tsi (timeserie
(range start (+ start window) (/ window intersteps))
0))
;value in[0 1] between each ref timestamps
;clockref (cf) ensure timestamps come from times
(defn intertimes[times ref]
(/
(cf (- (cf (time-as-value times)) (time-as-value ref)))
(- (backward (time-as-value ref)) (time-as-value ref))))
;linera interpolation
;clockref (cf) ensure timestamps come from times
;slice erase last uncomputable point
(defn linear-interpolation[times ref]
(def
dt (intertimes times ref)
p1 ref
p2 (backward ref))
(slice (+
p1
(cf (* (- p2 p1) (cf dt)))) [0 -1]))
;cubic interpolation
;clockref (cf) ensure timestamps come from times
;slice erase last uncomputable point
(defn cubic-interpolation[times ref]
(def
dt (intertimes times ref)
p0 (forward ref)
p1 ref
p2 (backward ref)
p3 (backward ref 2))
(slice (+
p1
(cf (* (+ (* -0.5 p0) (* 0.5 p2)) (cf dt)))
(* (+ p0 (* -2.5 p1) (* 2 p2) (* -0.5 p3)) (sqr dt))
(* (+ (* -0.5 p0) (* 1.5 p1) (* -1.5 p2) (* 0.5 p3)) (pow dt 3)))
[0 -1]))
;show reference and interpolations
[
(label "cos(14πt²)+t²" ts)
(label "linear" (linear-interpolation tsi ts))
(label "cubic" (cubic-interpolation tsi ts))
]
Simple Linear
If all timestamps of reference timeserie are included into requested timestamps we can use this simplifed version:
(defn linear-interpolation[times ref]
(+
ref
(*
(- (backward ref) ref)
(/
(- (time-as-value times)) (time-as-value ref)))
(- (backward (time-as-value ref)) (time-as-value ref))))
VWAP
Compute Volum Weighted Average Price since beginning of the day by creating vwap function.
Graph versus trade prices.
Use graphsample to get same graph but with less data transmission.
(defn vwap [source code date]
(label (str "vwap " code)
(/
(+
(*
(timeserie @trade-price source code date)
(timeserie @trade-volume source code date)))
(+
(timeserie @trade-volume source code date)))))
[
(graphsample 1920 (timeserie @trade-price "bitstamp" "BTC" 2018-02-17))
(graphsample 1920 (vwap "bitstamp" "BTC" 2018-02-17))
]
TWAP
Compute Time-Weighted Average Price on a 10 minutes sliding window by creating twap function.
Use graphsample to get same graphical representation but with less data transmission.
;time weighted over w
(defn twap[w ts]
(label (str "twap" w)
(/
(sliding w + (* (forward ts) (delta (time-as-value ts))))
(sliding w + (delta (time-as-value ts))))))
[
(graphsample 1920 (timeserie @trade-price "bitstamp" "BTC" 2018-02-16))
(graphsample 1920 (twap
10m
(timeserie @trade-price "bitstamp" "BTC" 2018-02-16)))
(graphsample 1920 (twap
1h
(timeserie @trade-price "bitstamp" "BTC" 2018-02-16)))
]
SMA EMA MACD
How to compute simple moving average (SMA), exponential moving average (EMA) and moving average convergence divergence (MACD) for any desired periods. (see Investopedia )
Use subsampled input timeserie to avoid too many points.
In this example we add 10K to MACD 12 26 9 to be graphed closer to other curves.
(defn sma[ts n]
"simple moving average for period n"
(/
(sliding n + ts)
n))
(defn ema[ts n]
"exponential moving average for period n"
(rts
(merge
(first (sma ts n))
(slice ts [(- n 1)])
)
(/ 2 (+ n 1))
(- 1 (/ 2 (+ n 1)))))
(defn macd[ts n1 n2 n3]
"MACD n1 n2 n3"
(ema
(-
(ema ts n1)
(ema ts n2)
)
n3))
(defn short[ts w]
"get last value for every w period"
(subsample w close ts))
;short name for serie
(def btc
(timeserie @trade-price "bitstamp" "BTC" 2018-02-16))
(def
w 10m ;on point per 10minutes
;define periods
n1 12
n2 26
n3 9)
; show curves
[
[
(subsample w open btc)
(subsample w high btc)
(subsample w low btc)
(subsample w close btc)
]
(label (str "sma " n1); nice label
(sma (short btc w) n1))
(label (str "ema " n1); nice label
(ema (short btc w) n1))
(label (str "ema " n2); nice label
(ema (short btc w) n2))
(label (str "macd " n1 " " n2 " " n3); nice label
(+ (macd (short btc w) n1 n2 n3) 10000))
]
Hayashi Yoshida
Compute Hayashi Yoshida correlation estimator over a slinding window.
Here is for example how Hayashi Yoshida could be implemented directly in LispTick. To have a nicer input we create a midserie function computing mid from bid and ask, keeping values in right time range.
We also compare result to trade price correlation which seems to be way more signifiant.
Use graphsample to get same graphical representation but with less data transmission.
(def
start 2020-03-02T09:10
stop 2020-03-02T17:30)
(defn midserie [code start stop]
(prune
(* 0.5
(+
(timeserie @bid-price "refinitiv" code start stop)
(timeserie @ask-price "refinitiv" code start stop)))))
(defn hayashi-yoshida [w x y]
(/
(sliding w +
(*
(delta x)
(delta y)))
(sqrt
(*
(sliding w + (* (delta x) (delta x)))
(sliding w + (* (delta y) (delta y)))))))
[
(label "1h mid correl"
(graphsample 1920
(hayashi-yoshida
1h
(midserie "BNPP.PA" start stop)
(midserie "SOGN.PA" start stop))))
(label "1h trade correl"
(graphsample 1920
(hayashi-yoshida
1h
(prune (timeserie @trade-price "refinitiv" "BNPP.PA" start stop))
(prune (timeserie @trade-price "refinitiv" "SOGN.PA" start stop)))))
]
Trade signing (aka trade-typing)
Trade signing determines if a trade is market buy or market sell. In other words whether a given trade is a liquidity-demander “buy” or liquidity-demander “sell”. Here we are presenting and implementing three popular conventions.
Tick Test
The tick test infers the direction of a trade by comparing its price to the price of the pre- ceding trade(s). The test classifies each trade into four categories: an uptick, a downtick, a zero-uptick, and a zero-downtick. A trade is an uptick (downtick) if the price is higher (lower) than the price of the previous trade. When the price is the same as the previous trade (a zero tick), if the last price change was a downtick, then the trade is a zero- downtick. A trade is classified as buy, 1 for our time series, if it occurs on an uptick or zero-uptick; otherwise it is classified as a sell, -1 for our time series.
Lee and Ready
Lee-Ready (see pdf) procedure, is based on both trade prices and quotes. First it uses what we call a “Side Test”, if trade price is higher than (lower than) best ask (best bid) it is a buy (a sell). For trade with price between bid and ask it uses Tick Test result.
Ellis, Michaely and O’Hara
Ellis, Michaely and O’Hara (see references) procedure, is also based on both trade prices and quotes. First it uses what we call a “Perfect Side Test”, if trade price is equal to best ask (best bid) it is a buy (a sell). For trade with price not perfectly equal to bid or ask it uses Tick Test result.
Implementation
Bitstamp is a European based cryptocurrency marketplace. Its web-socket API allows to retrieve lots of information for each trade and quote rarely available on standard market place. For example, for each trade we have timestamp, in second, price, amount… But we also have id, buy-order-id, sell-order-id and more important for this example: type, which is buy or sell. So having the real direction of the trade will allow us to compute an accuracy score for each method.
In this example we implement each convention and compare their accuracy each hour since beginning of the day.
(defn tick-test[source code start stop]
(*
(clockref
(one
(timeserie @trade-price source code start stop))
0 ;sign when no delta price available
)
(keep
(sign
(delta
(timeserie @trade-price source code (- start 1D) stop)))
not= 0)))
(defn lee-ready[source code start stop]
"Lee and Ready"
(merge
(keep
(side-test source code start stop)
not= 0)
(+
(clockref
(keep
(side-test source code start stop)
= 0))
(tick-test source code start stop))))
(defn emo[source code start stop]
"Ellis, Michaely, and O'Hara"
(merge
(keep ;use perfect-side if sign different from 0
(perfect-side-test source code start stop)
not= 0)
(+ ;use tick-test when perfect-side sign is 0
(clockref
(keep
(perfect-side-test source code start stop)
= 0))
(tick-test source code start stop))))
(defn perfect-side-test[source code start stop]
"-1 if trade = bid, 1 if trade = ask, 0 for all others"
(+ ;sign from ask plus sign from bid
(+ 1 (* -1 (abs (sign ; 1 if trade = ask 0 else
(-
(clockref ;only at trade time
(timeserie @trade-price source code start stop)
0 ;sign when no ask available
)
(timeserie @ask-price source code start stop))))))
(+ -1 (abs (sign ; -1 if trade = bid 0 else
(-
(clockref ;only at trade time
(timeserie @trade-price source code start stop)
0 ;sign when no bid available
)
(timeserie @bid-price source code start stop)))))))
(defn ask-side[source code start stop]
(sign
(+
1 ; only prices < ask will have 0 sign, others 1
(sign
(-
(clockref ;only at trade time
(timeserie @trade-price source code start stop)
0 ;sign when no ask available
)
(timeserie @ask-price source code start stop))))))
(defn bid-side[source code start stop]
(sign
(+
-1 ; only prices > bid will have 0 sign, others -1
(sign
(-
(clockref ;only at trade time
(timeserie @trade-price source code start stop)
0 ;sign when no bid available
)
(timeserie @bid-price source code start stop))))))
(defn side-test[source code start stop]
(+
(ask-side source code start stop)
(bid-side source code start stop)))
(def
source "bitstamp"
code "BTC"
d1 2019-02-21
d2 (+ d1 1D))
(def ref
(* -1 (timeserie @"sign" source code d1 d2)))
[
(label "EMO"
(/
(+ (time-truncate 1h (subsample 1h count
(keep (- ref (emo source code d1 d2)) = 0))))
(+ (time-truncate 1h (subsample 1h count ref)))))
(label "Lee-Ready"
(/
(+ (time-truncate 1h (subsample 1h count
(keep (- ref (lee-ready source code d1 d2)) = 0))))
(+ (time-truncate 1h (subsample 1h count ref)))))
(label "TickTest"
(/
(+ (time-truncate 1h (subsample 1h count
(keep (- ref (tick-test source code d1 d2)) = 0))))
(+ (time-truncate 1h (subsample 1h count ref)))))
(label "[0%-100%] scale" (timeserie d1 1 d1 0 d2 0))
]
Hampel Filter
The Hampel filter is a member of the class of decision filters that replaces the central value in the data window with the median if it lies far enough from the median to be deemed an outlier. A more in depth presentation can be found in Generalized Hampel Filters.
; Hampel Filter on timeserie ts
; w is half window size
; n number of sigmas to check for outliers
(defn hampel-filter[ts w n]
; filter point by point ab is [value histogram]
; median & sd done directly on histogram (last ab)
(defn one-hampel[ab]
(def x0 (median (last ab)))
; scale factor for Gaussian distribution
(def s0 (* 1.4826 (sd (last ab))))
(cond
(> (abs (- (first ab) x0)) (* n s0)) x0
(first ab)))
(map one-hampel
(sync ;synchronize value at t and histo at t+w
ts
(clockref
; 2*w+1 sliding window histogram
(sliding (+ (* 2 w) 1) hist
(backward ts w))))))
; Apply on real data
(def start 2018-06-06)
(def stop 2018-06-24)
(def temperature
(timeserie @"Temp_int" "mellisphera" "Ruche000" start stop))
[
(label "T°raw" temperature)
(label "T°filtered"
(hampel-filter temperature 5 1.2))
]
Finance algorithms applied to hive monitoring
Applying finance time series algorithms to hive temperature allows honey bee health monitoring.
When honey bees are sick, they are unable to keep inside hive temperature constant. So if inside temperature is strongly correlated to outside temperature honey bees are probably sick.
Here we are using Hayashi-Yoshida to compute correlation between in and out T° as they are asynchronous time series.
Then we smooth result using time-weighted average. If result is bigger than a defined threshold we switch alert ON.
Data gracefully provided by Mellisphera
(defn hayashi-yoshida [w x y]
"correlation over w"
(/
(sliding w +
(*
(delta x)
(delta y)))
(sqrt
(*
(sliding w + (* (delta x) (delta x)))
(sliding w + (* (delta y) (delta y)))))))
(defn twap[w ts]
"time weighted over w"
(/
(sliding w + (* (forward ts) (delta (time-as-value ts))))
(sliding w + (delta (time-as-value ts)))))
(defn alert[in out whayashi wtwap th]
"alert when weighted average correlation >threshold"
(sign (- (twap wtwap (hayashi-yoshida whayashi in out)) th)))
;parameters definition
(def start 2018-05-15)
(def stop 2018-07-07)
(def out
(timeserie @"Temp_ext" "mellisphera" "Ruche000" start stop))
(def in
(timeserie @"Temp_int" "mellisphera" "Ruche000" start stop))
(defn norm[ts]
(* 5 (+ 1 ts)))
(def w1 12h)
(def w2 6h)
(def threshold 0.83)
[
(label "T°hive" in)
(label "Ok" (norm (* -1 (alert in out w1 w2 threshold))))
(label "T°out" out)
(label "Alert" (norm (alert in out w1 w2 threshold)))
(label "Raw correlation" (hayashi-yoshida w1 in out))
(label "Smoothed correlation" (twap w2 (hayashi-yoshida w1 in out)))
]
Meteorology Statistics
In this example, we compute autocorrelation on some meteorology data. Once again, we use Hayashi-Yoshida between a timeserie and delayed copys of the same timeserie to get a full picture of autocorrelation.
Data comes from MeteoNet, it spans 3 years (2016-2018) with approximatly one point ever 6 minutes. We compute autocorrelation for 3 values:
- Humidity (humidité in french)
- Temperature (température)
- Pressure (pression)
WARNING can be slow to graph (>20s): 181 Millions points explored.
- each serie is 250.000 points long
- each step needs full serie and its delayed copy (x2)
- there are 121 steps for each serie (step is 24 minutes)
- there are 3 series
So a total of 250.000x2x121x3 = 181 Millions points to explore.
What can we conclude?
- Maximum value at 0s because every timeserie is perfectly correlated with itself.
- Each curve is perfectly sysmetrical because stop is the opposite of start.
- Humidity and temperature show same behavior with maximum correlation at ±24h and mininum at ±12h.
- This is due to night and day cycle, current temperature (or humidity) variation (increase or decrease) has the same direction as 24h ago.
- Current temperature (or humidity) variation has the opposite direction as 12h ago.
- Pressure autocorrelation is more “funny” with 2 maxima at ±24h and ±12h, this is a nice evidence of atmospheric tide.
(defn hayashi-yoshida [x y]
"full serie Hayashi Yoshida correlation estimator"
(/
(vget (last (+
(*
(delta (prune x))
(delta (prune y))))))
(sqrt
(vget (last (*
(last (+ (sqr (delta (prune x)))))
(last (+ (sqr (delta (prune y)))))
))))))
(defn autocorrelation[ts w]
"Autocorrelation with a w lag"
(cond ;trick to hit cache
(< w 0D) (hayashi-yoshida (- ts w) ts)
(hayashi-yoshida ts (+ ts w))))
;used arguments
(def
start 2016-01-01
stop 2018-12-31
biard "86027001"; aéroport Biard
station biard
step 24m
)
(defn ts[field]
"MeteoNet station field full timeserie"
(timeserie field "meteonet" station start stop))
(defn tscor[field step]
"timeserie autocorrelation from -24 to +24h step by step"
(label (str field)
(map-reduce-arg
(fn[w] (autocorrelation (ts field) w))
all
(range -24h 24h step))))
;show autocorrel for humidity, temperature and pressure
[
(label "humidité" (tscor @"hu" step))
(label "température" (tscor @"t" step))
(label "pression" (tscor @"psl" step))
]