Upgrade to Pro — share decks privately, control downloads, hide ads and more …

SOSTAT 2021: Time Series Analysis

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
Tweet

More Decks by Dr. Abbie Stevens

Other Decks in Science

Transcript

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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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”

    View full-size slide

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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  27. Wavelets are fantastic for gravitational waves!
    IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 28
    See also: “Q-transform”

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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 affects spectra
    IAA-SOSTAT 2021 ☆ Abbie Stevens, MSU & UMich 32

    View full-size slide

  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

    View full-size slide

  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)

    View full-size slide

  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

    View full-size slide

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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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 (figure); 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

    View full-size slide

  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

    View full-size slide

  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:
    signal adds,
    noise cancels
    CXY
    (ν)=X(ν)Y*(ν)
    real
    imaginary
    43
    real
    imaginary

    View full-size slide

  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

    View full-size slide

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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

  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)

    View full-size slide

  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

    View full-size slide

  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

    View full-size slide

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

    View full-size slide