Dr. Abbie Stevens
December 01, 2021
210

# SOSTAT 2021: Time Series Analysis

Introduction lecture to time series analysis in astronomy for the 2nd Severo Ochoa School on Statistics, Data Mining, and Machine Learning in Granada, Spain: https://www.granadacongresos.com/sostat2021

Tutorial in Jupyter notebooks: https://github.com/abigailStev/timeseries-tutorial

More on Dr. Abbie Stevens: https://abigailstevens.com/

## Dr. Abbie Stevens

December 01, 2021

## Transcript

1. Time series analysis
Dr. Abbie Stevens
Michigan State University & University of Michigan
[email protected]
@abigailStev
github.com/abigailstev

2. Outline
•Intro to the Fourier domain
•Evenly-spaced and irregular time series
•Power spectra/periodograms
•Coffee break and tutorial time
•Wavelets, Hilbert-Huang transform
•Spectrograms/dynamical power spectra
•Examples of astronomical signals
•More coffee and more tutorial time
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 2

3. The Fourier series
Any function can be represented as a sum of sines and
cosines (with some coefficients, which may also be functions)
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 3
f(x) = ∑ An
cos(n 𝜋 x) + ∑ Bn
sin(n 𝜋 x)
n=0 n=1
∞ ∞
Slide adapted from J. McIver

4. The Fourier series
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 4
Image credit:
J. Belk via
Wikimedia
Example: first four Fourier
approximation terms for a
square wave
The more terms you add, the
closer the approximation gets

5. The Fourier transform
Take a periodic or well-bounded function (of time or
space) by projecting f(t) onto an orthogonal basis of
sines and cosines
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 5
F(𝜈) = ∫ f(t) e-i2𝜋𝜈t dt
f(t) ⇾ F(𝜈) [or f(𝜈)]
^
f(t) = ∫ F(𝜈) ei2𝜋𝜈t d𝜈
Fourier transform
Inverse Fourier
transform
Slide adapted from J. McIver

6. The Fourier transform
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich GIF: L. Vieira via Wikimedia 6
Slide adapted from J. McIver
f(t) ⇾ F(𝜈) [or f(𝜈)]
^
Think of it like
decomposing the time
series function into its
component frequencies

7. The Fourier transform
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
Problem
solution
solve
(hard)
Transformed
problem
Transformed
solution
solve
(easy)
8
Fourier
transform
inverse Fourier
transform
The Fourier transform and inverse Fourier transform
make “Fourier pairs”

8. Sampling effect: Nyquist frequency
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 9
The more terms you add, the
closer the approximation gets

9. Sampling effect: Nyquist frequency
•In practice, cannot add terms forever!
•Highest frequency you can sample:
𝜈Nyquist
= 1/(2*dt) = d𝜈/2
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 10
The more terms you add, the
closer the approximation gets

10. Positive and negative Fourier frequencies
• Mirrored about 0 Hz
• In an array, often arranged: (0, pos., Nyquist, neg.)
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 11
Image: J. VanderPlas 2018

11. Evenly-sampled time series
• Signal period << observation length
• In X-rays and gamma-rays, we count photons. It’s very possible
to have zero counts -- “sparse” light curves are common
• Instead of saving light curves with lots of zeroes, we use
event lists – can select a dt multiple of the detector’s dt
• In optical, bright-enough sources mean you detect flux above
background in every time bin (e.g., every 30 seconds)
• Kepler, K2, TESS can give evenly sampled time series
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
1700 1702 1704 1706 1708 1710
2000 4000 6000 8000 10
4
1.2×10
4
Count/sec
Time (s)
Start Time 12339 7:28:14:566 Stop Time 12339 7:29:32:683
Bin time: 0.7812Eï02 s
12

12. Applying Fourier transforms to data
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
1016 1018 1020 1022 1024
5000 104 1.5×104
Count/sec
Time (s)
Start Time 10168 18:16:52:570 Stop Time 10168 18:17:08:180
Bin time: 0.1562Eï01 s
Time domain
Light curve
Frequency/Fourier domain
Power density spectrum
FOURIER
TRANSFORM2
Light curve broken into equal-length chunks or segments,
take power spectrum of each chunk, average those together
x(t)→X(ν)
P(ν)=X(ν)X*(ν)
=|X(ν)|2
13

13. Irregularly sampled time series
• Many (most?) astronomical time series will be irregularly sampled
Ø Signal period > observation length
Ø Weather
Ø Visibility
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 14
Reference: J. VanderPlas 2018, “Understanding the Lomb-Scargle Periodogram”

14. Lomb-Scargle periodogram saves the day!
• Multiplies your Fourier signal X(ν) with the Fourier transform
of the sampling window ⇾ “convolution”
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 15
Reference: J. VanderPlas 2018, “Understanding the Lomb-Scargle Periodogram”
When interpreting:
beware many noisy
peaks and harmonics
due to the sampling
window convolved
with noise in the data

15. X-ray variability: Hard to see by eye
1014 1016 1018 1020 1022
5000 104 1.5×104
Count/sec
Time (s)
CYGNUS_Xï1
Start Time 10168 18:16:52:578 Stop Time 10168 18:17:02:547
Bin time: 0.3125Eï01 s
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
1696 1698 1700 1702 1704
4000 5000 6000 7000
Count/sec
Time (s)
GRS1915+105
Start Time 12339 7:28:14:582 Stop Time 12339 7:28:24:542
Bin time: 0.4000Eï01 s
Light
curves
Power
density
spectra
Noise: Cygnus X-1 Signal: GRS 1915+105
16

16. QPOs → Damped harmonic oscillators
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
y = cos(⍵t)
17
I will discuss a little QPO physics
in the next lecture later today

17. QPOs → Damped harmonic oscillators
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
y = cos(⍵t) x e-bt
b=0
b=0.08
18

18. IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
y = cos(⍵t) x e-bt
b=0
b=0.08
b=0.22
19
QPOs → Damped harmonic oscillators

19. QPOs → Damped harmonic oscillators
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
y = cos(⍵t) x e-bt
b=0
b=0.08
b=0.22
b=0.5
20

20. QPOs → Damped harmonic oscillators
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
y = cos(⍵t) x e-bt
b=0
b=0.08
b=0.22
b=0.5
b=1.0
The stronger the damping,
the wider the peak
21

21. Poisson noise (“white noise”)
Poisson noise from
counting photons;
power-law slope=0
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 22

22. Poisson noise + Lorentzian QPO
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 23

23. Sampling effect: Aliasing
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 24
Resonance
between signal
freq. and
sample freq.
gives false
feature
Video credit: Honda
Windowing: similar
false feature due to
segment (“window”)
length

24. Tutorial time
•In the SOSTAT2021 Jupyter hub:
tutorials/abigail_timeseries/timeseries_workbook.ipynb
•On my GitHub:
https://github.com/abigailStev/timeseries-tutorial
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 25

25. Wavelets
• Fourier products don’t have
an intrinsic way to tell time
resolution (i.e., when in the
light curve the signal is
present)
• Wavelets easily represent a
signal in the time domain
and in the frequency
domain
Resource: “A really friendly guide to wavelets”, C. Valens, 1999
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 26

26. Wavelets aren’t great for everything
•Averaged power spectra (~50+ segments) follow a
chi-squared distribution with 2 degrees of freedom,
about the underlying true power spectrum
• Errors are statistically well-defined and well-understood (and
easy to compute!)
•Wavelets do not follow such a well-defined and well-
known distribution
•No clear, easy way to assess statistical significance of a
signal (which is one of the things we often want to do)
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 27

27. Wavelets are fantastic for gravitational waves!
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 28

28. Time (s)
40 41 42 43 44 45 46 47 48 49 50
Frequency (Hz)
1
2
3
4
5
6
7
8
9
10
Gaussian smoothing amplitude
0
5
10
15
20
25
30
Hilbert-Huang transform
• Frequency-domain product designed for data that are non-
stationary and non-linear
• Like an instantaneous Fourier transform → gives (some) time
localization!
• Error from standard
deviation of (1000+)
simulations
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich See Su+15 & refs therein for application to black hole data
0.01 0.1 1 10 100
10 100
Frequency(Hz)
Power (Leahy)
Hilbert spectrum from Su+15 of a 4Hz QPO
29

29. Elapsed time (in 64 s segments)
Power (rms2/Hz)
Elapsed time (in 64 s segments; not continuous)
Spectrogram (dynamical power spectrum)
•Evolution of a power spectrum in time
•Instead of averaging together, plot in colormap
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 30

30. Elapsed time (in 64 s segments)
Power (rms2/Hz)
Elapsed time (in 64 s segments; not continuous)
Good Time Intervals (GTIs)
• Jumps in the spectrogram below are cut out
• Looks like tophat windows in big light curve, but the segments we use for
the periodogram and averaging are much smaller than the window length
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 31

31. Things to sometimes worry about
• Deadtime occurs with X-ray detectors if your (bright) source is
emitting photons faster than your detector can handle. Once some
chip of the detector has detected a photon, it cannot detect another
photon until it reads out its existing photon detection through the
electronics.
ØMeasurable as deviation from expected Poisson noise power-
law at high frequencies in the power spectrum
ØAccumulates over an observation; for a few ks observation,
could have several seconds of deadtime to adjust the exposure
time by
• Pile-up is a sibling of deadtime that aﬀects spectra
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 32

32. Things to sometimes worry about
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 33
Slide adapted from D. Huppenkothen; credit: NuSTAR Observatory Guide
‘Deadtime’ ⇾ the detector is effectively dead for the brief readout period.
NuSTAR

33. Things to sometimes worry about
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 34
Slide adapted from D. Huppenkothen; credit: Huppenkothen & Bachetti 2021 (in press)
Huppenkothen & Bachetti (under review)

34. Poisson noise (a form of “white noise”)
• In power spectra: P~ν𝛼, 𝛼=0
• Not correlated on any timescale
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 35

35. Flicker noise
• In power spectra: P~ν𝛼, 𝛼=-1
• Correlated on medium-short timescales (short “memory”)
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 36

36. Red noise (random walk)
• In power spectra: P~ν𝛼, 𝛼=-2
• Correlated on long timescales (long “memory”)
• Ornstein-Uhlenbeck: ~red noise + friction: tends
towards a mean value over long time
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 37

37. Beware of red noise!
• Cannot apply standard peak-finding algorithms, since those
assume white noise (see Vaughan & Uttley ‘06)
• Bigger issue for SMBHs
than stellar BHs
due to timescales
White/Poisson noise
Red noise
Smith+18b
Fourier frequency
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 38

38. Red noise vs signals
−2 0 2 4 6 8 10 12
time (yr)
V mag

● ●

15.2 15.0 14.8
(a)
−2 0 2 4 6 8 10 12
time (yr)
V mag

● ●

15.2 15.0 14.8
(b)

(c)
PG 1302-102, CRTS data
Vaughan+16 (figure); Liu+18
Data looks periodic!
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 39

39. Red noise vs signals
−2 0 2 4 6 8 10 12
time (yr)
V mag

● ●

15.2 15.0 14.8
(a)
−2 0 2 4 6 8 10 12
time (yr)
V mag

● ●

15.2 15.0 14.8
(b)

(c)
PG 1302-102, CRTS data
Vaughan+16 (figure); Liu+18
Uneven sampling, gappy
data, only ~1.5 cycles
Data looks periodic!
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 40

40. Red noise vs signals
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
−2 0 2 4 6 8 10 12
time (yr)
V mag

● ●

15.2 15.0 14.8
(a)
−2 0 2 4 6 8 10 12
time (yr)
V mag

● ●

15.2 15.0 14.8
(b)

(c)
Vaughan+16 (ﬁgure); Liu+18
PG 1302-102, CRTS data
Also including LINEAR data
(but it isn’t)
Sampling a random red
noise process in same
way can look like a
“periodic” signal
Uneven sampling, gappy
data, only ~1.5 cycles
Data looks periodic!
41

41. Red noise vs signals
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
−2 0 2 4 6 8 10 12
time (yr)
V mag

● ●

15.2 15.0 14.8
(a)
−2 0 2 4 6 8 10 12
time (yr)
V mag

● ●

15.2 15.0 14.8
(b)

(c)
Vaughan+16 (figure); Liu+18
Data looks periodic!
Uneven sampling, gappy
data, only ~1.5 cycles
PG 1302-102, CRTS data
Sampling a random red
noise process in same
way can look like a
“periodic” signal
Also including LINEAR data
When in doubt, simulate!
Also: claimed periodicity in
J0045+41 disproven by
Barth & Stern ‘18
42

42. Cross spectra
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
x(t)→X(ν) for a narrow energy band
y(t)→Y(ν) for a broad-energy reference band
As you average
segments together:
noise cancels
CXY
(ν)=X(ν)Y*(ν)
real
imaginary
43
real
imaginary

43. Cross spectra
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
x(t)→X(ν) for a narrow energy band
y(t)→Y(ν) for a broad-energy reference band
Cospectrum: real part of the cross spectrum (see Bachetti+15,
Bachetti+Huppenkothen 18 and Huppenkothen+Bachetti 18 for statistical details)
Note: for X(ν)=Y(ν),
cospectrum = cross amplitude
= power spectrum
44
Also used:
amplitude of the
cross spectrum

44. What other science cases might use
these timing analysis tools?
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich

45. Strohmayer 2001
0.1 1 10 100
2
4
6
8
Type-A
0.1 1 10 100
101
Leahy Power
Type-B
0.1 1 10 100
Frequency [Hz]
101
Type-C
Stevens
Mo1a+17a
QPOs in black holes and neutron stars
• High-frequency: 100’s Hz
• Hot Keplerian blobs in inner disk?
• Low-frequency: ~0.01-10’s Hz
• Precession of corona/hot flow?
Magnetic warps in disk?
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 46

46. • Short-timescale
variability changes
on long timescales
(spectral state-
dependent)
• Short-timescale
variability is energy-
dependent
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
Corona-dominated state
Mostly-corona-dominated state
Disk-dominated state
Few months
for full
spectral state
transition
Mostly-disk-dominated state
If you want to read more, see the power colours paper by Heil+15a 47
QPOs in black holes and neutron stars

47. • Short-timescale
variability changes
on long timescales
(spectral state-
dependent)
• Short-timescale
variability is energy-
dependent
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
A B
C D
2-60 keV
6.5-13.1 keV 13.1-60 keV
2-6.5
keV
Homan+01 48
QPOs in black holes and neutron stars

48. Pulsations in neutron stars
• Spin-down: decreasing spin frequency
(e.g., losing rotational energy to the
environment)
• Spin-up: increasing spin frequency (e.g.,
accreting material and thus increasing
angular momentum)
• Glitch: sudden change in spin frequency
(due to superfluid NS core?)
• Seen in residuals of frequency or pulse
timing
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
-400
-300
-200
-100
0
Timing residuals (ms)
(a)
-100 -50 0 50 100
Days from MJD = 53067.1
-3.750
-3.745
-3.740
-3.735
ν (10-10 Hz s-1)

(d)
-100 -50 0 50 100
0
1
2
3
4
5
6
Δν (μHz)
(c)
-100 -50 0 50 100
-200
-100
0
100
200
Timing residuals (ms)
(b)
Espinoza+11 49

49. Stellar pulsations
• Cepheid variables, RR Lyrae stars, Delta Scuti variables,
Blahzko effect (long-period
modulation of the periodicity)
• Period-luminosity relation
makes them standard candles
used as “cosmic distance ladder”
• Slow enough (periods > hours)
that time-domain photometry
is often used
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich Image via APOD, credit: R. Vanderbei, ESA/Gaia/DPAC 50

50. Asteroseismology (“starquakes”)
• Understanding the internal structure of stars using their
brightness oscillations
• Convective zone excites oscillations
• Fourier analysis of light curves:
often see many harmonics
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich Info thanks to online slides by T. Bedding and refs therein
Aerts+19
51

51. QPOs in active galactic nuclei (AGN)
• ~1 hr “periodicity”, 91ks observation
• RE J1034+396 is a narrow-line Seyfert 1
AGN
• Saw 16 ‘cycles’ (periods) in one
uninterrupted observation!
• Evenly-sampled time bins
• Signal attributed to high-freq. QPO
• If at innermost stable circular orbit,
MBH
~7x106-1x107 M☉
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich
Gierlinski+08
Alston+14
52

52. Smith+18b
• 44 day low-freq. QPO in KIC 9650712
• NLS1 in original Kepler field
• 30-minute cadence over 3.5 years: ~30 cycles
• Tested periodicity via simulations (Uttley+02) and Lomb-
Scargle periodogram
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 53
QPOs in active galactic nuclei (AGN)

53. • Open-source timing and spectral-timing software
(Astropy affiliated package)
Øhttps://docs.stingray.science
• Stingray: Python library of analysis tools
• HENDRICS: shell scripting interface
• DAVE: graphical user interface
• Tutorials in Jupyter notebooks
• Huppenkothen, Bachetti, Stevens et al. 2019, ApJ & JOSS
• Google Summer of Code students in 2016-2021
Stingray
Please remember to cite
software in your papers!
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 54

54. Time series resources
Evenly-sampled time series:
• Uttley et al. 2014, “X-ray Reverberation Around Accreting
Black Holes”
Irregularly-sampled time series:
• VanderPlas 2018, “Understanding the Lomb-Scargle
Periodogram”
Wavelets:
• C. Valens 1999, “A really friendly guide to wavelets”
Software tools:
• Stingray, Lightkurve, GWpy, astropy.timeseries
IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 55

55. IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 56