<jupyter_start><jupyter_text>Lecture 14: Introduction to Time Series Analysis  This notebook was developed by [<NAME>](http://faculty.washington.edu/ivezic/) for the 2021 data science class at the University of Sao Paulo and it is available from [github](https://github.com/ivezic/SaoPaulo2021/blob/main/notebooks/Lecture14.ipynb).Note: this notebook contains code developed by <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME> and many others.  Modeling Toolkit For Time Series Analysis  This notebook includes:[Fourier Analysis](fourier)[Discrete Fourier Transform](dft) [An example: estimating period for variable stars](periods)[Digital Filtering](filtering) [Bonus topic: Analysis of Stochastic Processes](stochastic)  IntroductionIn this lecture, we are going to learn a few fundamental concepts and tools for analyzing time series data. From an astronomical perspective, this field has predominantly focused on unevenly sampled data, low signal-to-noise data, and heteroscedastic errors. The most common applications range from the detection of variability and periodicity to the treatment of nonperiodicvariability and searches for localized events.  Despite being set in the context of time series, many tools and results are readily applicable in other domains, which we will see later in this notebook. We will extend this toolkit by introducing tools for analysisin the frequency domain, such as digital filtering.  The main data sets used in this chapter include light curves obtained by the [LINEAR survey](https://arxiv.org/abs/1505.02082) (Lincoln Near-Earth Asteroid Research).  Main concepts for time series analysisThe time series discussed here will be limited to two-dimensional scalar data sets: pairs of random variables, ($t_1$, $y_1$),...,($t_N$; $y_N$), with no assumptions about the sampling of the time coordinate t. In many ways, analysis methods discussed here are closely related to the parameter estimation and model selection problems discussed in the context of regression. However, unlike regression problems where different y measurements are typically treated as independent randomvariables, **for time-series models the value of $y_{i+1}$ directly depends on the preceding value $y_i$**.    The main tasks of time series analysis are:  1. To characterize the presumed temporal correlation between different values of y, including its significance.  2. To forecast (predict) future values of y.  Import Functions and Packages<jupyter_code>import numpy as np
from matplotlib import pyplot as plt

from astroML.datasets import fetch_rrlyrae_templates

from scipy.signal import fftconvolve

from scipy import fftpack
from matplotlib import mlab

from astroML.datasets import fetch_LIGO_large

from scipy.fftpack import fft
from scipy.stats import norm

from astroML.fourier import PSD_continuous
from astroML.plotting import setup_text_plots
setup_text_plots(usetex=True)<jupyter_output><empty_output><jupyter_text>1. Fourier analysis [Go to top](toc)Fourier analysis plays a major role in the analysis of time series data. In Fourier analysis, generalfunctions are represented or approximated by integrals or sums of simpler trigonometric functions.   For periodic functions, such as periodic light curvesin astronomy, it is often true that a relatively small number of terms (less than 10) suffices to reachan approximation precision level similar to the measurement precision.  Some useful conceptsThe **Fourier transform** of function h(t) is defined as$$H(f)=\int_{-\infty}^{\infty} h(t) exp(-i2\pi ft)dt$$with **inverse transformation**$$h(t)=\int_{-\infty}^{\infty} H(t) exp(-i2\pi ft)df$$where t is time and f is frequency (for time in seconds, the unit for frequency is hertz, or Hz; theunits for H(f) are the product of the units for h(t) and inverse hertz.  For a real function h(t), H(f) is in general a complex function.  In the special case when h(t) is an even function such that h(-t) = h(t), H(f) is real and even as well. For example, theFourier transform of a pdf of a zero-mean Gaussian $N(0, \sigma)$ in the time domain is a Gaussian$H(f) = exp(-2\pi^2 \sigma^2 f^2)$ in the frequency domain.  When the time axis of an arbitrary functionh(t) is shifted by $\Delta t$, then the Fourier transform of $h(t + \Delta t)$ is$$\int_{-\infty}^{\infty}h(t+\Delta t)exp(-i2\pi ft)dt = H(f)exp(i2\pi f\Delta t)$$Therefore, the **Fourier transform of a Gaussian $N(\mu, \sigma)$** is$$H_{Gauss}(f)=exp(-2\pi^2\sigma^2f^2)[cos(2\pi f\mu)+isin(2\pi f \mu)]$$This is known as **"white noise"** since there is no frequencydependence (also known as "thermal noise" or Johnson's noise).  An important quantity in time series analysis is the one-sided **power spectral density (PSD)** function(or power spectrum) defined for $0 \leq f < \infty$ as$$PSD(f) \equiv |H(f)|^2+|H(-f)|^2$$The PSD gives the amount of power contained in the frequency interval between f and f + df.  The total power is the same whether computed in the frequency or the time domain:$$P_{tot} \equiv \int_0^{\infty}PSD(f)df=\int_{-\infty}^{\infty}|h(t)|^2dt$$This result is known as **Parseval's theorem**.  1.1 Fourier Reconstruction of RR-Lyrae TemplatesBelow is an example of a truncated Fourier representation of an RR Lyrae light curve. As we can see, the more terms that are included in the sum, the better is the resulting approximation.  Import a RR Lyrae templateFrom astroML.datasets, we take the RR Lyrae (variable star) light curve as an example. RR Lyrae has a periodical pulsation.<jupyter_code>templates = fetch_rrlyrae_templates()
x, y = templates['115r'].T<jupyter_output><empty_output><jupyter_text>Show Fourier reconstruction resultsThe thick dashed line shows the true curve; the gray lines show the approximation based on 1, 3, and 8 Fourier modes (sinusoids).<jupyter_code>fig = plt.figure(figsize=(5, 5))
fig.subplots_adjust(hspace=0)

kvals = [1, 3, 8]
subplots = [311, 312, 313]

for (k, subplot) in zip(kvals, subplots):
    ax = fig.add_subplot(subplot)

    # Use FFT to fit a truncated Fourier series
    y_fft = np.fft.fft(y)
    y_fft[k + 1:-k] = 0
    y_fit = np.fft.ifft(y_fft).real

    # plot the true value and the k-term reconstruction
    ax.plot(np.concatenate([x, 1 + x]),
            np.concatenate([y, y]), '--k', lw=2)
    ax.plot(np.concatenate([x, 1 + x]),
            np.concatenate([y_fit, y_fit]), color='gray')

    label = "%i mode" % k
    if k > 1:
        label += 's'

    ax.text(0.02, 0.1, label, ha='left', va='bottom',
            transform=ax.transAxes)

    if subplot == subplots[-1]:
        ax.set_xlabel('phase')
    else:
        ax.xaxis.set_major_formatter(plt.NullFormatter())

    if subplot == subplots[1]:
        ax.set_ylabel('amplitude')
    ax.yaxis.set_major_formatter(plt.NullFormatter())

    ax.set_xlim(0, 2)
    ax.set_ylim(1.1, -0.1)

plt.show()<jupyter_output><empty_output><jupyter_text>1.2 Convolution TheoremA convolution of two functions $a(t)$ and $b(t)$ is given by$$(a*b)(t) \equiv \int_{-\infty}^{\infty} a(t')b(t-t')dt'$$In the above equation, the function a canbe thought of as the "convolving pattern" of the measuring apparatus, and the function $b$ is thesignal. In practice, we measure the convolved (or smoothed) version of our signal, $[a * b](t)$, andseek to uncover the original signal $b$ using the presumably known $a$.  The convolution theorem states that if $h = a * b$, then the Fourier transforms of $h$, $a$, and $b$ arerelated by their pointwise products:$$H(f)=A(f)B(f)$$Thus a convolution of two functions is transformed into a simple multiplication of the associatedFourier representations. Therefore, to obtain $b$, we can simply take the inverse Fourier transformof the ratio $H(f)=A(f)$. In the absence of noise, this operation is exact.    For further explanation, we will see a schematic representation of the convolution theorem below.  2. Discrete Fourier Transform[Go to top](toc)In practice, data are always discretely sampled. When the spacing of the time interval is constant,the discrete Fourier transform is a powerful tool.Although in astronomy, temporal data are rarely sampled with uniform spacing (LIGO data, however, as a counterexample, will be discussed in the later notebook), uniformly sampled data is a good place to start because we have fast algorithm for this situation.  When computing the Fourier transform for discretely and uniformly sampled data, the Fourierintegrals from $$H(f)=\int_{-\infty}^{\infty} h(t) exp(-i2\pi ft)dt$$$$h(t)=\int_{-\infty}^{\infty} H(t) exp(-i2\pi ft)df$$are translated to sums. Here we have t for time and f for frequency.  The discrete Fourier transform of the vector of values hj is a complex vector of length N definedby$$H_k=\sum_{j=0}^{N-1}h_j exp[-i2\pi jk/N]$$where k=0,...,(N-1).  The corresponding inverse discrete Fourier transform is defined by$$h_j=\sum_{j=0}^{N-1}H_k exp[-i2\pi jk/N]$$where j=0,...,(N-1). Unlike the continuous transforms, here the units for $H_k$ are the same asthe units for $h_j$ . Given $H_k$, we can represent the function described by $h_j$ as a sum of sinusoids,as was done in the first figure.  2.1 The Nyquist sampling theoremLet us define $h(t)$ to be band limited if $H(f) = 0$ for $|f| > f_c$, where $f_c$ is the band limit, or theNyquist critical frequency. If $h(t)$ is band limited, then there is some "resolution" limit in $t$ space,$t_c = 1/(2f_c)$ below which h(t) appears "smooth." When $h(t)$ is band limited, then according to the Nyquist sampling theorem we can exactly reconstruct h(t) from evenly sampled data when $\Delta t < t_c$,as$$h(t)=\frac{\Delta t}{t_c} \sum_{k=-\infty}^{k=\infty} h_k \frac{sin[2\pi f_c (t-k\Delta t)]}{2\pi f_c (t-k\Delta t)}$$This result is known as the Whittaker-Shannon, or often just Shannon, interpolation formula(or "sinc-shifting" formula).  Note that the summation goes to infinity, but also that the termmultiplying $h_k$ vanishes for large values of $|t - k\Delta t|$.  For example $h(t) = sin(2 \pi t/P )$ has a periodP and is band limited with $f_c = 1/P$. If it is sampled with $\Delta t$ not larger than $P/2$, it can be fullyreconstructed at any $t$ (it is important to note that this entire discussion assumes that there is nonoise associated with sampled values $h_j$ ). On the other hand, when the sampled function h(t) is notband limited, or when the sampling rate is not sufficient (i.e., $\Delta t > t_c$), an effect called **"aliasing"**prevents us from exactly reconstructing h(t). In such a case, all of the powerspectral density from frequencies $|f| > f_c$ is aliased (falsely transferred) into the $-f_c < f < f_c$range. The aliasing can be thought of as inability to resolve details in a time series at a finer detailthan that set by $f_c$. The figure below will show us the effect of aliasing.  2.2 The Fast Fourier Transform (FFT)The **Fast Fourier transform (FFT)** is an algorithm for computing discrete Fourier transforms inO(N logN) time, rather than O(N2) using a naive implementation. The algorithmic details forthe FFT can be found in NumRec. The speed of FFT makes it a widespread tool in the analysisof evenly sampled, high signal-to-noise ratio, time series data.  An example below shows such analysis for a function with a single dominant frequency: a sine wave whose amplitude is modulated by a Gaussian.  Generate the data<jupyter_code>np.random.seed(1)

tj = np.linspace(-25, 25, 512)
hj = np.sin(tj)
hj *= norm(0, 10).pdf(tj)<jupyter_output><empty_output><jupyter_text>Show result with data and calculated PSDThe discrete Fourier transform (bottom panel) for two noisy data sets shown in the top panel. For 512 evenly sampled times t (dt = 0.977), points are drawn from h(t) = a + sin(t)G(t), where G(t) is a Gaussian N(mu = 0,sigma = 10). Gaussian noise with sigma = 0.05 (top data set) and 0.005 (bottom data set) is added to signal h(t). The value of the offset a is 0.15 and 0, respectively. The discrete Fourier transform is computed as described in this section. For both noise realizations, the correct frequency f = (2pi)-1 ~ 0.159 is easily discernible in the bottom panel. Note that the height of peaks is the same for both noise realizations. The large value of abs(H(f = 0)) for data with larger noise is due to the vertical offset.<jupyter_code># plot the results
fig = plt.figure(figsize=(5, 3.75))
fig.subplots_adjust(hspace=0.35)
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)

offsets = (0, 0.15)
colors = ('black', 'gray')
linewidths = (1, 2)
errors = (0.005, 0.05)

for (offset, color, error, linewidth) in zip(offsets, colors,
                                             errors, linewidths):
    # compute the PSD
    err = np.random.normal(0, error, size=hj.shape)
    hj_N = hj + err + offset
    fk, PSD = PSD_continuous(tj, hj_N)

    # plot the data and PSD
    ax1.scatter(tj, hj_N, s=4, c=color, lw=0)
    ax1.plot(tj, 0 * tj + offset, '-', c=color, lw=1)
    ax2.plot(fk, PSD, '-', c=color, lw=linewidth)

# vertical line marking the expected peak location
ax2.plot([0.5 / np.pi, 0.5 / np.pi], [-0.1, 1], ':k', lw=1)

ax1.set_xlim(-25, 25)
ax1.set_ylim(-0.1, 0.3001)

ax1.set_xlabel('$t$')
ax1.set_ylabel('$h(t)$')

ax1.yaxis.set_major_locator(plt.MultipleLocator(0.1))

ax2.set_xlim(0, 0.8)
ax2.set_ylim(-0.101, 0.801)

ax2.set_xlabel('$f$')
ax2.set_ylabel('$PSD(f)$')

plt.show()<jupyter_output><empty_output><jupyter_text>3. An example: estimating period for variable stars[Go to top](toc)The code below is essentially the code that generates Figure 10.17 in the textbook, seehttps://www.astroml.org/book_figures/chapter10/fig_LINEAR_LS.htmlThat code:* downloads a dataset of variable stars collected with the [LINEAR](https://en.wikipedia.org/wiki/Lincoln_Near-Earth_Asteroid_Research) telescope,* selects six of them (as given by the `ids` variable),* computes the Lomb-Scargle periodogram in the `compute_best_frequencies` function,* finds the best frequency by looking for the highest peak in the periodiogram, and* folds the light curves and plots the results.We will below modify that code to do the following:1. Find the frequencies for the first 100 stars in the LINEAR datasets (i.e., those with IDs in `data.ids[:100]`, where `data` is the variable containing the LINEAR dataset as defined below).2. Compute the periods corresponing to the frequencies above, and plot the histogram of periods using 50 bins in range $0 < T < 25$ hours.  3. Generate two 6-panel figures analogous to Figure 10.17, for 6 stars (arbitrarily chosen) from the highest and second highest peak, respectivelly.<jupyter_code># Figure 10.17 code

# Author: <NAME>
# License: BSD
#   The figure produced by this code is published in the textbook
#   "Statistics, Data Mining, and Machine Learning in Astronomy" (2013)
#   For more information, see http://astroML.github.com
#   To report a bug or issue, use the following forum:
#    https://groups.google.com/forum/#!forum/astroml-general
import numpy as np
from matplotlib import pyplot as plt

from astroML.decorators import pickle_results
from astroML.time_series import search_frequencies, lomb_scargle, MultiTermFit
from astroML.datasets import fetch_LINEAR_sample

#----------------------------------------------------------------------
# This function adjusts matplotlib settings for a uniform feel in the textbook.
# Note that with usetex=True, fonts are rendered with LaTeX.  This may
# result in an error if LaTeX is not installed on your system.  In that case,
# you can set usetex to False.
from astroML.plotting import setup_text_plots
setup_text_plots(fontsize=8, usetex=True)

#------------------------------------------------------------
# Compute the best frequencies
@pickle_results('LINEAR_LS.pkl')
def compute_best_frequencies(ids, n_eval=10000, n_retry=5, generalized=True):
    results = {}
    for i in ids:
        t, y, dy = data[i].T
        print(" - computing power for %i (%i points)" % (i, len(t)))
        kwargs = dict(generalized=generalized)
        omega, power = search_frequencies(t, y, dy, n_eval=n_eval,
                                          n_retry=n_retry,
                                          LS_kwargs=kwargs)
        results[i] = [omega, power]

    return results

# plot phased light curves 
def plot_lightcurves(results, ids):
    #------------------------------------------------------------
    # Plot the phased light-curves
    fig = plt.figure(figsize=(5, 6.5))
    fig.subplots_adjust(hspace=0.1, bottom=0.06, top=0.94, left=0.12, right=0.94)

    for i in range(len(ids)):
        # get the data and best-fit angular frequency
        t, y, dy = data[ids[i]].T
        omega, power = results[ids[i]]
        omega_best = omega[np.argmax(power)]
        print(" - omega_0 = %.10g" % omega_best)

        # do a fit to the first 4 Fourier components
        mtf = MultiTermFit(omega_best, 4)
        mtf.fit(t, y, dy)
        phase_fit, y_fit, phased_t = mtf.predict(1000, return_phased_times=True)

        # plot the phased data and best-fit curves
        ax = fig.add_subplot(321 + i)
        ax.errorbar(phased_t, y, dy, fmt='.k', ecolor='gray',
                    lw=1, ms=4, capsize=1.5)
        ax.plot(phase_fit, y_fit, '-b', lw=2)

        ax.set_xlim(0, 1)
        ax.set_ylim(plt.ylim()[::-1])
        ax.yaxis.set_major_locator(plt.MaxNLocator(4))

        ax.text(0.03, 0.04, "ID = %i" % ids[i], ha='left', va='bottom',
                transform=ax.transAxes)
        ax.text(0.03, 0.96, "P = %.2f hr" % (2 * np.pi / omega_best * 24.),
                ha='left', va='top',
                transform=ax.transAxes)

        ylim = ax.get_ylim()
        ax.set_ylim(ylim[0], ylim[0] + 1.1 * (ylim[1] - ylim[0]))

        if i < 4:
            ax.xaxis.set_major_formatter(plt.NullFormatter())
        if i % 2 == 0:
            ax.set_ylabel('mag')

        if i in (4, 5):
            ax.set_xlabel('phase')

    return fig 
# Load the dataset
data = fetch_LINEAR_sample()
# select a subset of the full dataset
ids = [14752041, 1009459, 10022663, 10025796, 11375941, 18525697]

# Compute the best frequencies
results = compute_best_frequencies(ids, n_eval=10000, n_retry=5)

# Plot the phased light-curves
myprettyplot = plot_lightcurves(results, ids)
# read data and process the first 100 stars
data = fetch_LINEAR_sample()
Nlightcurve = 100

# if you are brave you can do all 7009! 
#  it took about 5 hours on my laptop
# if you want to do that, uncomment this line: 
# Nlightcurve = data.ids[:-1].size

ids = data.ids[:Nlightcurve]
results = compute_best_frequencies(ids, n_eval=10000, n_retry=5)

omegas = np.zeros(len(ids))
for i in range(len(ids)):
    # get the data and best-fit angular frequency
    t, y, dy = data[ids[i]].T
    omega, power = results[ids[i]]
    omega_best = omega[np.argmax(power)]
    omegas[i] = omega_best
    print(" - Star #%4d: omega_0, T_0 = %.10g, %.10g" % (i, omega_best, (2 * np.pi / omega_best * 24.)))

# Convert frequencies to periods (in hours):
TT = (2 * np.pi / omegas * 24.)
# Construct and plot the histogram of periods:
n, bins, patches = plt.hist(TT, 50, range=(0,25))
plt.xlabel('Period (hour)')
plt.ylabel('N')
plt.show()
# Now pick 6 ids close to the first and second peaks (determined by eyeballing of the above diagram):
ids_peak1 = ids[(TT >= 3) & (TT < 3.5)][:6]
ids_peak1
ids_peak2 = ids[(TT >= 14) & (TT < 16)][:6]
ids_peak2
# And plot the two diagrams:
p1 = plot_lightcurves(results, ids_peak1)
p2 = plot_lightcurves(results, ids_peak2)<jupyter_output>- omega_0 = 10.20062676
 - omega_0 = 10.29120901
 - omega_0 = 10.0363784
 - omega_0 = 9.613121674
 - omega_0 = 10.00077894
 - omega_0 = 9.448396134<jupyter_text>4. Digital Filtering [Go to top](toc)  Introduction**Digital filtering** aims to reduce noise in time series data, or to compress data. Common examplesinclude low-pass filtering, where high frequencies are suppressed, high-pass filtering, where lowfrequencies are suppressed, passband filtering, where only a finite range of frequencies is admitted,and a notch filter, where a finite range of frequencies is blocked. We will use a few examples to illustrate the most commonapplications of filtering.  Fourier analysis is one of the most useful tools for performing filtering. Numerous other techniques can be found in signal processing literature, including approaches based on the wavelets discussed in the modeling toolkit notebook.     We emphasize that filtering always decreases the information content of data (despite makingit appear less noisy). As we have already learned throughout previous chapters, when modelparameters are estimated from data, raw (unfiltered) data should be used. In some sense, this isan analogous situation to binning data to produce a histogram-while very useful for visualization, estimates of model parameters can become biased if one is not careful.  Import packages and dataIn this notebook, we are going to explore the astroML.filters. We mainly use the Savitzky-Golay and the Wiener filters. The spectrum of a white dwarf data imported for this notebook is from Sloan Digital SkySurvey (SDSS).<jupyter_code>import numpy as np
from matplotlib import pyplot as plt

from astroML.fourier import PSD_continuous
from astroML.datasets import fetch_sdss_spectrum

from scipy import optimize, fftpack, interpolate
from scipy.signal import savgol_filter
from astroML.fourier import IFT_continuous
from astroML.filters import wiener_filter
from astroML.filters import min_component_filter<jupyter_output><empty_output><jupyter_text>4.1 Low-pass filtersThe power spectrum for common Gaussian noise is at and will extend to frequencies as high asthe Nyquist limit, $f_N = 1=(2\Delta t)$. If the data are band limited to a lower frequency, $f_c < f_N$, thenthey can be smoothed without much impact by suppressing frequencies $|f| > f_c$. Given a filterin frequency space, $\Phi (f)$, we can obtain a smoothed version of data by taking the inverse Fouriertransform of$$\hat{Y}(f)=Y(f) \Phi(f)$$where $Y(f)$ is the discrete Fourier transform of data. At least in principle, we could simply set$\Phi(f)$ to zero for $|f| > f_c$, but this approach would result in ringing (i.e., unwanted oscillations)in the signal. Instead, the optimal filter for this purpose is constructed by minimizing the MISEbetween $\hat{Y}(f)$ and $Y(f)$ (for detailed derivation see NumRec) and is called the **Wiener filter**:$$\Phi(f) = \frac{P_S(f)}{P_S(f)+P_N(f)}$$Here $P_S(f)$ and $P_N(f)$ represent components of a two-component (signal and noise) fit to thePSD of input data, $PSD_Y (f) = P_S(f) + P_N(f)$, which holds as long as the signal and noise areuncorrelated. We will see how the filtering works in the example below.  Create the noisy dataWe are going to generate a set of noisy data as the input signal, on which we apply filters. The figure below shows the input data (200 evenly spaced points) with a narrow Gaussian peak centered at x = 20.<jupyter_code># Create the noisy data
np.random.seed(5)
N = 2000
dt = 0.05

t = dt * np.arange(N)
h = np.exp(-0.5 * ((t - 20.) / 1.0) ** 2)
hN = h + np.random.normal(0, 0.5, size=h.shape)

# Plot the results
N = len(t)
Df = 1. / N / (t[1] - t[0])
f = fftpack.ifftshift(Df * (np.arange(N) - N / 2))
HN = fftpack.fft(hN)

fig = plt.figure(figsize=(5, 3.75))
fig.subplots_adjust(wspace=0.05, hspace=0.35,
                    bottom=0.1, top=0.95,
                    left=0.12, right=0.95)

# First plot: noisy signal
ax = fig.add_subplot(111)
ax.plot(t, hN, '-', c='gray')
ax.plot(t, np.zeros_like(t), ':k')
ax.text(0.98, 0.95, "Input Signal", ha='right', va='top',
        transform=ax.transAxes, bbox=dict(fc='w', ec='none'))

ax.set_xlim(0, 90)
ax.set_ylim(-0.5, 1.5)

ax.xaxis.set_major_locator(plt.MultipleLocator(20))
ax.set_xlabel(r'$\lambda$')
ax.set_ylabel('flux')<jupyter_output><empty_output><jupyter_text>Set up the Wiener filterWe fit a model to the PSD consisting of the sum of a gaussian and white noise using **Wiener filter**. We will see this method in the later sections in this notebook.<jupyter_code>N = 2000

Df = 1. / N / dt
f = fftpack.ifftshift(Df * (np.arange(N) - N / 2))
HN = fftpack.fft(hN)

# apply the Wiener filter
h_smooth, PSD, P_S, P_N, Phi = wiener_filter(t, hN, return_PSDs=True)<jupyter_output><empty_output><jupyter_text>Set up the Savitzky-Golay filterWe use a fourth-order **Savitzky-Golay** filter with a window size of $\Delta \lambda = 10$ to filter the vales. The Savitzky-Golay filter is a very simple but powerful method as a low-pass filter. It fits low-order polynomialsto data (in the time domain) using sliding windows (it is also known as the least-squares filter).For a detailed discussion, see NumRec.<jupyter_code># apply the Savitzky-Golay filter
h_sg = savgol_filter(hN, window_length=201, polyorder=4, mode='mirror')<jupyter_output><empty_output><jupyter_text>Show filtered signalPlot below shows noisy signal after filtering. * Result from Wiener filter is shown in black.* Result from Savitzky-Golay filter is shown in gray.  The Gaussian peak at x=20 is clearly seen in both curves.<jupyter_code># Second plot: filtered signal
ax = plt.subplot(111)
ax.plot(t, np.zeros_like(t), ':k', lw=1)
ax.plot(t, h_smooth, '-k', lw=1.5, label='Wiener')
ax.plot(t, h_sg, '-', c='gray', lw=1, label='Savitzky-Golay')

ax.text(0.98, 0.95, "Filtered Signal", ha='right', va='top',
        transform=ax.transAxes)
ax.legend(loc='upper right', bbox_to_anchor=(0.98, 0.9), frameon=False)

ax.set_xlim(0, 90)
ax.set_ylim(-0.5, 1.5)

ax.xaxis.set_major_locator(plt.MultipleLocator(20))
ax.set_xlabel(r'$\lambda$')
ax.set_ylabel('flux')<jupyter_output><empty_output><jupyter_text>Show filtered PSD* The upper panel shows the input power spectral density (PSD) distribution.  * The lower panel shows the Wiener-filtered power spectral density (PSD) distributions.  The two curves in the upper panel represent two-component fit to PSD given by equation$$\Phi(f) = \frac{P_S(f)}{P_S(f)+P_N(f)}$$<jupyter_code># Plot the results
N = len(t)
Df = 1. / N / (t[1] - t[0])
f = fftpack.ifftshift(Df * (np.arange(N) - N / 2))
HN = fftpack.fft(hN)

fig = plt.figure(figsize=(5, 3.75))
fig.subplots_adjust(wspace=0.05, hspace=0.35,
                    bottom=0.1, top=0.95,
                    left=0.12, right=0.95)

# Third plot: Input PSD
ax = fig.add_subplot(211)
ax.scatter(f[:N // 2], PSD[:N // 2], s=9, c='k', lw=0)
ax.plot(f[:N // 2], P_S[:N // 2], '-k')
ax.plot(f[:N // 2], P_N[:N // 2], '-k')

ax.text(0.98, 0.95, "Input PSD", ha='right', va='top',
        transform=ax.transAxes)

ax.set_ylim(-100, 3500)
ax.set_xlim(0, 0.9)

ax.yaxis.set_major_locator(plt.MultipleLocator(1000))
ax.xaxis.set_major_locator(plt.MultipleLocator(0.2))
ax.set_xlabel('$f$')
ax.set_ylabel('$PSD(f)$')

# Fourth plot: Filtered PSD
ax = fig.add_subplot(212)
filtered_PSD = (Phi * abs(HN)) ** 2
ax.scatter(f[:N // 2], filtered_PSD[:N // 2], s=9, c='k', lw=0)

ax.text(0.98, 0.95, "Filtered PSD", ha='right', va='top',
        transform=ax.transAxes)

ax.set_ylim(-100, 3500)
ax.set_xlim(0, 0.9)

ax.yaxis.set_major_locator(plt.MultipleLocator(1000))
ax.xaxis.set_major_locator(plt.MultipleLocator(0.2))
ax.set_xlabel('$f$')
ax.set_ylabel('$PSD(f)$')<jupyter_output><empty_output><jupyter_text>Wiener Filter and kernel smoothing ConnectionThere is an interesting connection between the kernel density estimation method discussed in the KDE notebookand Wiener filtering. By the convolution theorem, the Wiener-filtered result is equivalent to theconvolution of the unfiltered signal with the inverse Fourier transform of $\Phi(f)$.  This convolution is equivalent to kernel density estimation. When Wiener filtering is viewed in this way, it effectively says that we believe the signal is as wide as the centralpeak, and the statistics of the noise are such that the minor peaks in thewings work to cancel out noise in the major peak.  Hence, the modeling of the PSD in the frequencydomain via $$\Phi(f) = \frac{P_S(f)}{P_S(f)+P_N(f)}$$ corresponds to choosing the optimal kernel width. Just as detailed modelingof the Wiener filter is not of paramount importance, the choice of kernel is not either.We will use the same data as the previous Wiener filter figure as an example to explore this connection.  Find effective kernel<jupyter_code># inverse fourier transform Phi to find the effective kernel
t_plot, kernel = IFT_continuous(f, Phi)<jupyter_output><empty_output><jupyter_text>perform kernel smoothingThis is faster in frequencyspace (i.e. using the standard Wiener filter above) but we will do it in the slow & simple way here to demonstrate the equivalence explicitly.<jupyter_code>kernel_func = interpolate.interp1d(t_plot, kernel.real)

t_eval = np.linspace(0, 90, 1000)
t_KDE = t_eval[:, np.newaxis] - t
t_KDE[t_KDE < t_plot[0]] = t_plot[0]
t_KDE[t_KDE > t_plot[-1]] = t_plot[-1]
F = kernel_func(t_KDE)

h_smooth = np.dot(F, hN) / np.sum(F, 1)<jupyter_output><empty_output><jupyter_text>Show kernel and smoothing results* The left panel shows the inverse Fourier transform of the Wiener filter $\Phi(f)$ applied in the sample data we used previously.  * The right panel shows the data smoothed by this kernel, which is equivalent to the Wiener filter smoothing in the previous figure.  By the convolution theorem, the Wiener-filtered result is equivalent to the convolution of the unfiltered signal withthe kernel shown above, and thus Wiener filtering and kernel density estimation (KDE) are directly related.<jupyter_code># Plot the results
fig = plt.figure(figsize=(10, 4.4))
fig.subplots_adjust(left=0.1, right=0.95, wspace=0.25,
                    bottom=0.15, top=0.9)

# First plot: the equivalent Kernel to the WF
ax = fig.add_subplot(121)
ax.plot(t_plot, kernel.real, '-k')
ax.text(0.95, 0.95, "Effective Wiener\nFilter Kernel",
        ha='right', va='top', transform=ax.transAxes)

ax.set_xlim(-10, 10)
ax.set_ylim(-0.05, 0.45)
ax.set_xlabel(r'$\lambda$')
ax.set_ylabel(r'$K(\lambda)$')

# Second axes: Kernel smoothed results
ax = fig.add_subplot(122)
ax.plot(t_eval, h_smooth, '-k', lw=1)
ax.plot(t_eval, 0 * t_eval, '-k', lw=1)
ax.text(0.95, 0.95, "Kernel smoothing\nresult",
        ha='right', va='top', transform=ax.transAxes)

ax.set_xlim(0, 90)
ax.set_ylim(-0.5, 1.5)

ax.set_xlabel('$\lambda$')
ax.set_ylabel('flux')

plt.show()<jupyter_output><empty_output><jupyter_text>4.2 High-pass filtersThe most common example of high-pass filtering in astronomy is baseline estimation in spectraldata. Unlike the case of low-pass filtering, here there is no universal filter recipe. Baseline estimationis usually the first step toward the estimation of model parameters (e.g. location, width, and strength of spectral lines). In such cases, the best approach might be full modeling andmarginalization of baseline parameters as nuisance parameters at the end of analysis.    A simple iterative technique for high-pass filtering, called **minimum component filtering**, is discussedin detail in WJ03. These are the **main steps**:1. Determine baseline: exclude or mask regions where signal is clearly evident and fit a baselinemodel (e.g., a low-order polynomial) to the unmasked regions.2. Get FT for the signal: after subtracting the baseline fit in the unmasked regions (i.e., a linearregression fit), apply the discrete Fourier transform.3. Filter the signal: remove high frequencies using a low-pass filter (e.g., Wiener filter), andinverse Fourier transform the result.4. Recombine the baseline and the filtered signal: add the baseline fit subtracted in step 2 tothe result from step 3. This is the minimum component filtering estimate of baseline.    In the next two examples, we will see the application of a minimum component filter to the spectrum of a white dwarf.  Example 1  Fetch the dataWe first fetch the spectrum data from SDSS database for use. The intermediate steps of the minimum component filter procedure applied to the spectrum of a white dwarf from the SDSS data set (mjd= 52199, plate=659, fiber=381).<jupyter_code># Fetch the spectrum from SDSS database
plate = 659
mjd = 52199
fiber = 381

data = fetch_sdss_spectrum(plate, mjd, fiber)

lam = data.wavelength()
spec = data.spectrum<jupyter_output><empty_output><jupyter_text>Pre-process the dataWavelengths we get are logorithmically spaced: we will work in log(lam).<jupyter_code>def preprocess(lam, spec):
    loglam = np.log10(lam)

    flag = (lam > 4000) & (lam < 5000)
    lam = lam[flag]
    loglam = loglam[flag]
    spec = spec[flag]

    lam = lam[:-1]
    loglam = loglam[:-1]
    spec = spec[:-1]
    
    return [lam, loglam, spec]

[lam, loglam, spec] = preprocess(lam, spec)<jupyter_output><empty_output><jupyter_text>Apply minimum component filtering steps First step: mask-out significant features<jupyter_code>feature_mask = (((lam > 4080) & (lam < 4130)) |
                ((lam > 4315) & (lam < 4370)) |
                ((lam > 4830) & (lam < 4900)))<jupyter_output><empty_output><jupyter_text>Second step: fit a line to the unmasked portion of the spectrum<jupyter_code>XX = loglam[:, None] ** np.arange(2)
beta = np.linalg.lstsq(XX[~feature_mask], spec[~feature_mask], rcond=None)[0]

spec_fit = np.dot(XX, beta)
spec_patched = spec - spec_fit
spec_patched[feature_mask] = 0<jupyter_output><empty_output><jupyter_text>Third step: Fourier transform the patched spectrum<jupyter_code>N = len(loglam)
df = 1. / N / (loglam[1] - loglam[0])
f = fftpack.ifftshift(df * (np.arange(N) - N / 2.))
spec_patched_FT = fftpack.fft(spec_patched)<jupyter_output><empty_output><jupyter_text>Fourth step: Low-pass filter on the transform<jupyter_code>filt = np.exp(- (0.01 * (abs(f) - 100.)) ** 2)
filt[abs(f) < 100] = 1

spec_filt_FT = spec_patched_FT * filt<jupyter_output><empty_output><jupyter_text>Fifth step: inverse Fourier transform, and add back the fit<jupyter_code>spec_filt = fftpack.ifft(spec_filt_FT)
spec_filt += spec_fit<jupyter_output><empty_output><jupyter_text>Show filter result* The top panel shows the input spectrum; the masked sections of the input spectrum are shown by thin lines (i.e., step 1 of the process). * The bottom panel shows the PSD of the masked spectrum, after the linear fit has been subtracted (gray line). * A simple low-pass filter (dashed line) is applied, and the resulting filtered spectrum (dark line) is used to construct the result shown in the next figure.<jupyter_code># plot results
fig = plt.figure(figsize=(5, 3.75))
fig.subplots_adjust(hspace=0.45)

ax = fig.add_subplot(211)
ax.plot(lam, spec, '-', c='gray')
ax.plot(lam, spec_patched + spec_fit, '-k')

ax.set_ylim(25, 110)
ax.set_xlim(4000, 5000)

ax.set_xlabel(r'$\lambda\ {\rm(\AA)}$')
ax.set_ylabel('flux')

ax = fig.add_subplot(212)
factor = 15 * (loglam[1] - loglam[0])
ax.plot(fftpack.fftshift(f),
        factor * fftpack.fftshift(abs(spec_patched_FT) ** 1),
        '-', c='gray', label='masked/shifted spectrum')
ax.plot(fftpack.fftshift(f),
        factor * fftpack.fftshift(abs(spec_filt_FT) ** 1),
        '-k', label='filtered spectrum')
ax.plot(fftpack.fftshift(f),
        fftpack.fftshift(filt), '--k', label='filter')

ax.set_xlim(0, 2000)
ax.set_ylim(0, 1.1)

ax.set_xlabel('$f$')
ax.set_ylabel('scaled $PSD(f)$')<jupyter_output><empty_output><jupyter_text>Example 2 (use package) Fetch the dataWe follow the same process to analyze the same white dwarf example. Here instead of explicitly applying steps, we call function *min_component_filter* to achieve the same goal.<jupyter_code>plate = 659
mjd = 52199
fiber = 381

data = fetch_sdss_spectrum(plate, mjd, fiber)

lam = data.wavelength()
spec = data.spectrum<jupyter_output><empty_output><jupyter_text>Pre-process the data<jupyter_code>[lam, loglam, spec] = preprocess(lam,spec)<jupyter_output><empty_output><jupyter_text>Apply minimum component filtering using functionThe function used here is *min_component_filter* in *astroML.filters* package.<jupyter_code>feature_mask = (((lam > 4080) & (lam < 4130)) |
                ((lam > 4315) & (lam < 4370)) |
                ((lam > 4830) & (lam < 4900)))

spec_filtered = min_component_filter(loglam, spec, feature_mask, fcut=100)<jupyter_output><empty_output><jupyter_text>Compute PSD of filtered and unfiltered versions<jupyter_code>f, spec_filt_PSD = PSD_continuous(loglam, spec_filtered)
f, spec_PSD = PSD_continuous(loglam, spec)<jupyter_output><empty_output><jupyter_text>Show filter result* The upper panel shows a portion of the input spectrum, along with the continuum computed via the minimum component filtering procedure described above (See the previous figure). * The lower panel shows the PSD for both the input spectrum and the filtered result.<jupyter_code># Plot the results
fig = plt.figure(figsize=(5, 3.75))
fig.subplots_adjust(hspace=0.45)

# Top panel: plot noisy and smoothed spectrum
ax = fig.add_subplot(211)
ax.plot(lam, spec, '-', c='gray', lw=1)
ax.plot(lam, spec_filtered, '-k')

ax.text(0.97, 0.93, "SDSS white dwarf\n %i-%i-%i" % (mjd, plate, fiber),
        ha='right', va='top', transform=ax.transAxes)

ax.set_ylim(25, 110)
ax.set_xlim(4000, 5000)

ax.set_xlabel(r'$\lambda\ {\rm (\AA)}$')
ax.set_ylabel('flux')

# Bottom panel: plot noisy and smoothed PSD
ax = fig.add_subplot(212, yscale='log')
ax.plot(f, spec_PSD, '-', c='gray', lw=1)
ax.plot(f, spec_filt_PSD, '-k')

ax.set_xlabel(r'$f$')
ax.set_ylabel('$PSD(f)$')
ax.set_xlim(0, 2000)<jupyter_output><empty_output><jupyter_text>Bonus topic: Analysis of Stochastic Processes[Go to top](toc)  IntroductionStochastic variability includes behavior that is not predictable forever as in the periodic case, butunlike temporally localized events, variability is always there. Typically, the underlying physics isso complex that we cannot deterministically predict future values.  Despite their seemingly irregular behavior, stochastic processes can be quantified. In this notebook, we will mainly discuss **autocorrelation**, **autoregressive** and **damped random walk** models in analyzing stochastic processes.  Import functionsWe will use functions from *astroML.time_series* and from *astroML.fourier*.<jupyter_code>from astroML.time_series import generate_power_law
from astroML.time_series import lomb_scargle, generate_damped_RW
from astroML.time_series import ACF_scargle, ACF_EK<jupyter_output><empty_output><jupyter_text>The autocorrelation and structure functionsOne of the main statistical tools for the analysis of stochastic variability is the **autocorrelationfunction**. It represents a specialized case of the correlation function of two functions, f(t) and g(t),scaled by their standard deviations, and defined at time lag $\Delta t$ as$$CF(\Delta t) = \frac{lim_{T \rightarrow \infty} \frac{1}{T} \int_{(T)}f(t)g(t+ \Delta t)dt}{\sigma_f \sigma_g}$$where $\sigma_f$ and $\sigma_g$ are standard deviations of f(t) and g(t), respectively. With this normalization, thecorrelation function is unity for $\Delta t$ = 0 (without normalization by standard deviation, the aboveexpression is equal to the covariance function). It is assumed that both f and g are statisticallyweakly stationary functions, which means that their mean and autocorrelation function (see below)do not depend on time. The correlation function yields information about the time delaybetween two processes. If one time series is produced from another one by simply shifting the timeaxis by tlag, their correlation function has a peak at $\Delta t = t_{lag}$.    With $f(t) = g(t) = y(t)$, the autocorrelation of y(t) defined at time lag $\Delta t$ is$$ACF(\Delta t) = \frac{lim_{T \rightarrow \infty} \frac{1}{T} \int_{(T)}y(t)y(t+ \Delta t)dt}{\sigma_y^2}$$The autocorrelation function yields information about the variable timescales present in a process. When y values are uncorrelated (e.g., due to white noise without any signal), $ACF(\Delta t) = 0$, except for $ACF(0) =1$. For processes that "retain memory" of previous states only for some characteristictime $\tau$, the autocorrelation function vanishes for $\Delta t \gg \tau$.     The autocorrelation function and the PSD of function y(t) (i.e. $PSD(f) \equiv |H(f)|^2 + |H(-f)|^2$) are Fourier pairs; thisfact is known as the Wiener-Khinchin theorem and applies to stationary random processes. Theformer represents an analysis method in the time domain, and the latter in the frequency domain.  The **structure function** is another quantity closely related to the autocorrelation function,$$SF(\Delta t) = SF_\infty [1-ACF(\Delta t)]^{\frac{1}{2}}$$where $SF_\infty$ is the standard deviation of the time series evaluated over an infinitely large timeinterval (or at least much longer than any characteristic timescale $\tau$).  Examples of stochastic processes: $1/f$ and $1/f^2$ processesFor a given autocorrelation function or PSD, the corresponding time series can be generated usingthe algorithm described in [On generating power law noise](http://adsabs.harvard.edu/full/1995A%26A...300..707T). Essentially, the amplitude of the Fourier transform is given by thePSD, and phases are assigned randomly; the inverse Fourier transform then generates time series.  The connection between the PSD and the appearance of time series of two power-law PSDs: $1/f$ and $1/f^2$ is illustrated in figure below.    1. Generate data for plottingLet us use 1024 data points, and let $\Delta t = 0.01$ in this case.<jupyter_code>N = 1024
dt = 0.01
factor = 100

t = dt * np.arange(N)
random_state = np.random.RandomState(1)<jupyter_output><empty_output><jupyter_text>2. Calculate PSD and compare resultsWe use *generate_power_law* in *astroML_time_series* package to apply the algorithm [On generating power law noise](http://adsabs.harvard.edu/full/1995A%26A...300..707T) on our generated data.  In the plot below shows examples of stochastic time series generated from power-law PSDs.* Plots of generated data (top panels) and the corresponding PSD (bottom panels) are shown.* Results from $1/f$ power-law (left panels) and from $1/f^2$ (right panel) are shown.* Dashed lines indicate input PSD, and solid lines are determined from time series shown in the top panels. The PSD normalization is such that both cases have similarpower at low frequencies. For this reason, the overall amplitudes (more precisely, the variance) ofthe two time series are similar. The power at high frequencies is much larger for the $1/f$ case, andthis is why the corresponding time series has the appearance of noisy data.  The structure function for the $1/f$ process is constant, and proportional to $t^{1/2}$ forthe $1/f^2$ process (remember that we defined structure function with a square root).<jupyter_code>fig = plt.figure(figsize=(10, 7.5))
fig.subplots_adjust(wspace=0.05)

for i, beta in enumerate([1.0, 2.0]):
    # Generate the light curve and compute the PSD
    x = factor * generate_power_law(N, dt, beta, random_state=random_state)
    f, PSD = PSD_continuous(t, x)

    # First axes: plot the time series
    ax1 = fig.add_subplot(221 + i)
    ax1.plot(t, x, '-k')

    ax1.text(0.95, 0.05, r"$P(f) \propto f^{-%i}$" % beta,
             ha='right', va='bottom', transform=ax1.transAxes)

    ax1.set_xlim(0, 10.24)
    ax1.set_ylim(-1.5, 1.5)

    ax1.set_xlabel(r'$t$')

    # Second axes: plot the PSD
    ax2 = fig.add_subplot(223 + i, xscale='log', yscale='log')
    ax2.plot(f, PSD, '-k')
    ax2.plot(f[1:], (factor * dt) ** 2 * (2 * np.pi * f[1:]) ** -beta, '--k')

    ax2.set_xlim(1E-1, 60)
    ax2.set_ylim(1E-6, 1E1)

    ax2.set_xlabel(r'$f$')

    if i == 1:
        ax1.yaxis.set_major_formatter(plt.NullFormatter())
        ax2.yaxis.set_major_formatter(plt.NullFormatter())
    else:
        ax1.set_ylabel(r'${\rm counts}$')
        ax2.set_ylabel(r'$PSD(f)$')<jupyter_output><empty_output><jupyter_text>Autocorrelation and structure function for evenly and unevenly sampled data Evenly sampled dataIn the case of **evenly sampled data**, with $t_i = (i-1)\Delta t$, the autocorrelation function of a discretely sampled $y(t)$ is defined as$$ACF(j) = \frac{\sum^{N-j}_{i=1} [(y_i-\bar{y})(y_{i+j}-\bar{y})]}{\sum^{N}_{i=1}(y_i-\bar{y})^2}$$The uncertainty of the mean is $$\sigma_{\bar{x}} = \frac{\sigma}{\sqrt{N}}[1+2 \sum^N_{j=1}(1-\frac{j}{N})ACF(j)]^{1/2}$$where $\sigma$ is the homoscedastic measurement error.   When a time series has a nonvanishing ACF, the uncertainty of its mean is larger than for an uncorrelated data set, where $\sigma_{\bar{x}} = \frac{s}{\sqrt{N}}$.  When data are unevenly sampled, the ACF cannot be computed using equations above.   Unevenly sampled dataFor the case of **unevenly sampled data**, <NAME> Krolik proposed the **"discrete correlation function" (DCF)** in an astronomical context (called the "slot autocorrelation function" in physics). See [The Discrete Correlation Function](https://ui.adsabs.harvard.edu/abs/1988ApJ...333..646E/abstract).  For discrete unevenly sampled data with homoscedastic errors, they defined a quantity$$UDCF_{ij} = \frac{(y_i-\bar{y})(g_j-\bar{g})}{[(\sigma_y^2-e_y^2)(\sigma_g^2-e_g^2)^{1/2}]}$$  where $e_y$ and $e_g$ are homoscedastic measurement errors for time series y and g. The associated timelag is $\Delta t_{ij} = ti - tj$. The discrete correlation function at time lag $\Delta t$ is then computed by binningand averaging $UDCF_{ij}$ over M pairs of points for which $\Delta t-\delta t=2 \leq \delta t_{ij} \leq \Delta t+\delta t/2$, where $\delta t$ isthe bin size. The bin size is a trade-of between accuracy of $DCF(\Delta t)$ and its resolution. Edelson and Krolik showed that even uncorrelated time series will produce values of the cross-correlation $DCF(\Delta t) \sim \pm 1/ \sqrt{M}$.  Scargle has developed different techniques to evaluate the discrete Fourier transform, correlationfunction and autocorrelation function of unevenly sampled time series (see [Studies in astronomical time series analysis. III.](https://www.osti.gov/biblio/5344858-studies-astronomical-time-series-analysis-iii-fourier-transforms-autocorrelation-functions-cross-correlation-functions-unevenly-spaced-data)).  Given an unevenly sampled time series, y(t), the essential steps of Scargle'sprocedure are as follows:1. Compute the generalized Lomb-Scargle periodogram for $y(t_i), i = 1,..., N$, namely $P_{LS}(\omega)$.2. Compute the sampling window function using the generalized Lomb-Scargle periodogram using $z(t_i) = 1, i = 1,..., N$, namely $P_{LS}^W(\omega)$.3. Compute inverse Fourier transforms for $P_{LS}(\omega)$ and $P^W_{LS}(\omega)$, namely $\rho(t)$ and $\rho^W(t)$, respectively.4. The autocorrelation function at lag t is $ACF(t) = \rho(t)/ \rho^W(t)$.  Edelson and Krolik's DCF method and the Scargle method DemonstrationWe will see an example of the use of Edelson and Krolik's DCF method and the Scargle method.  1. Generate time-series dataWe use do 1000 days worth of magnitudes for this sample data.<jupyter_code>t = np.arange(0, 1E3)
z = 2.0
tau = 300
tau_obs = tau / (1. + z)

np.random.seed(6)
y = generate_damped_RW(t, tau=tau, z=z, xmean=20)

# randomly sample 100 of these
ind = np.arange(len(t))
np.random.shuffle(ind)
ind = ind[:100]
ind.sort()
t = t[ind]
y = y[ind]

# add errors
dy = 0.1
y_obs = np.random.normal(y, dy)<jupyter_output><empty_output><jupyter_text>2. Compute ACF via Scargle methodWe calculate the ACF following Scargle method process.<jupyter_code>C_S, t_S = ACF_scargle(t, y_obs, dy,
                       n_omega=2. ** 12, omega_max=np.pi / 5.0)

ind = (t_S >= 0) & (t_S <= 500)
t_S = t_S[ind]
C_S = C_S[ind]<jupyter_output><empty_output><jupyter_text>3. Compute ACF via E-K methodWe calculate ACF using Edelson and Krolik's DCF method.<jupyter_code>C_EK, C_EK_err, bins = ACF_EK(t, y_obs, dy, bins=np.linspace(0, 500, 51))
t_EK = 0.5 * (bins[1:] + bins[:-1])<jupyter_output><empty_output><jupyter_text>4. Compare resultsExample of the autocorrelation function for a stochastic process is shown below. * The top panel shows a simulated light curve generated using a damped random walk model. * The bottom panel shows the corresponding autocorrelation function computed using Edelson and Krolik's DCF method and the Scargle method. * The solid line shows the input autocorrelation function used to generate the light curve.  The two methods produce similar results. Errors are easier to compute for the DCF method and this advantageis crucial when fitting models to the autocorrelation function.<jupyter_code># Plot the results
fig = plt.figure(figsize=(5, 5))

# plot the input data
ax = fig.add_subplot(211)
ax.errorbar(t, y_obs, dy, fmt='.k', lw=1)
ax.set_xlabel('t (days)')
ax.set_ylabel('observed flux')

# plot the ACF
ax = fig.add_subplot(212)
ax.plot(t_S, C_S, '-', c='gray', lw=1,
        label='Scargle')
ax.errorbar(t_EK, C_EK, C_EK_err, fmt='.k', lw=1,
            label='Edelson-Krolik')
ax.plot(t_S, np.exp(-abs(t_S) / tau_obs), '-k', label='True')
ax.legend(loc=3)

ax.plot(t_S, 0 * t_S, ':', lw=1, c='gray')

ax.set_xlim(0, 500)
ax.set_ylim(-1.0, 1.1)

ax.set_xlabel('t (days)')
ax.set_ylabel('ACF(t)')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11707018081520074, 'fraction_numerical': 0.025232403718459494, 'mean_word_length': 4.4036211084124535, 'pattern_counts': {'":': 0, '<': 169, '<?xml version=': 0, '>': 166, 'https://': 7, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '35', 'dir': 'jupyter-structured-clean-dedup', 'id': '154620', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 16911, 'n_tokens_neox': 15536, 'n_words': 5728, 'path': 'lectures/notes/Lecture5-time-series.ipynb', 'repo_name': 'uw-astro/astr-598a-win22', 'size': '48994'}