Time Series

Lecture 11

February 18, 2026

Review

GLMs

  • Generalized Linear Models: use linear models to represent alternative distributional assumptions.
  • \[y_i \sim \mathcal{D}(\theta), \quad f(\theta) = \sum_j \beta_j x_i^j\]
  • Focus on logistic and Poisson regressions.
  • Can be more challenging to interpret coefficients.

Time Series

Time Series and Environmental Data

Many environmental datasets involve time series: repeated observations over time:

\[X = \{X_t, X_{t+h}, X_{t+2h}, \ldots, X_{t+nh}\}\]

More often: ignore sampling time in notation,

\[X = \{X_1, X_2, X_3, \ldots, X_n\}\]

Figure 1: Example of time series; trapped lynx populations in Canada.

Spectral Analysis

Often time series have periodic components. Can recover these using regression: \[y_t = A \cos(\omega t) + B \sin(\omega t) + C + \varepsilon_t, \qquad \varepsilon \sim N(0, \sigma^2).\]

Tide Gauge Data

Code 
function load_data(fname)
    date_format = "yyyy-mm-dd HH:MM"
    # this uses the DataFramesMeta package -- it's pretty cool
    return @chain fname begin
        CSV.File(; dateformat=date_format)
        DataFrame
        DataFrames.rename(
            "Time (GMT)" => "time", "Predicted (m)" => "harmonic", "Verified (m)" => "gauge"
        )
        @transform :datetime = (Date.(:Date, "yyyy/mm/dd") + Time.(:time))
        select(:datetime, :gauge, :harmonic)
        @transform :weather = :gauge - :harmonic
        @transform :month = (month.(:datetime))
    end
end

dat = load_data("data/surge/norfolk-hourly-surge-2015.csv")

p1 = plot(dat.datetime, dat.gauge; ylabel="Gauge Measurement (m)", label="Observed", legend=:topleft, xlabel="Date/Time", bottom_margin=5mm, left_margin=5mm, right_margin=5mm, size=(1100, 500))

Figure 2: 2015 tide gauge data from the Norfolk, VA tide gauge.

Example: Monthly Periods

If we want to look at the tide gauge data (which has an hourly sampling period and 8760 samples) with a monthly period, then \[\begin{gathered} 2\pi / \omega = 8760 / 30 = 292 \\ \omega = 2\pi / 292 \end{gathered}\]

and we get: \([A, B, C, \sigma\)] = $[-0.012, -0.022, 0.516, 0.34].

How To Identify Interesting Periods

Spectral analysis: Based on Fourier decomposition to switch between time domain and frequency domain.

Estimating Fourier Coefficients

Let \(\omega_j = 2\pi j / T\). Then:

\[\begin{aligned} \hat{A} &= \frac{2}{T} \sum_t y_t \cos(\omega_j t) \\ \hat{B} &= \frac{2}{T} \sum_t y_t \sin(\omega_j t) \\ \hat{C} &= \bar{Y} = \frac{1}{T} \sum_t y_t. \end{aligned}\]

Periodogram

Can then define the periodogram \[I_T(\omega) = \frac{1}{2\pi T} \left|\sum_{t=1}^T y_t e^{i\omega t}\right|^2.\]

Larger values of \(I_T\) mean those frequencies contribute more to the overall signal.

Periodogram: Tide Gauge

Code 
T = nrow(dat) # get number of samples (this is hourly in this case), so 8760
ps = DSP.periodogram(dat.gauge; fs=T) # compute periodogram with T sampling frequency
plot(ps.freq, ps.power, xlabel="Period (Hours)", ylabel="Power Spectrum", legend=false, linewidth=3, yaxis=:log)
Figure 3: Periodogram for the NOrfolk Tide Gauge Data

The biggest spikes are at period 0 (no period; this is just a drift without harmonic behavior) and period 711 (roughly 30 days; aligns with the lunar cycle).

Sometimes periodograms are computed with respect to frequency, rather than period: recover period by multiplying by the sampling frequency.

Autocorrelation

Temporal Dependence

Dependence: History or sequencing of the data matters

\[p(y_t) = f(y_1, \ldots, y_{t-1})\]

Serial dependence captured by autocorrelation:

\[\varsigma(i) = \rho(y_t, y_{t-i}) = \frac{\text{Cov}[y_t, y_{t-i}]}{\mathbb{V}[y_t]} \]

Figure 4: Autocorrelation for the Lynx data.

Stationarity

I.I.D. Data: All data drawn from the same distribution, \(y_i \sim \mathcal{D}\).

Equivalent for time series \(\mathbf{y} = \{y_t\}\) is stationarity.

  • Strict stationarity: \(\{y_t, \ldots, y_{t+k-1}\} \sim \mathcal{D}\) for all \(t\).
  • Weak stationarity: \(\mathbb{E}[X_1] = \mathbb{E}[X_t]\) and \(\text{Cov}(X_1, X_k) = \text{Cov}(X_t, X_{t+k-1})\) for all \(t, k\).

Stationary vs. Non-Stationary Series

(a) Stationary vs. Non-stationary time series
(b)
Figure 5

Autocorrelation

(a) Lag-1 autocorrelation for an AR model and Gaussian noise.
(b)
(c)
(d)
Figure 6

Autocorrelation and Linear Regression

Code 
# tds_riverflow_loglik: function to compute the log-likelihood for the tds model
# θ: vector of model parameters (coefficients β₀ and β₁ and stdev σ)
# tds, flow: vectors of data
tds = let
    fname = "data/tds/cuyaTDS.csv" # CHANGE THIS!
    tds = DataFrame(CSV.File(fname))
    tds[!, [:date, :discharge_cms, :tds_mgL]]
end

function tds_riverflow_loglik(θ, tds, flow)
    β₀, β₁, σ = θ # unpack parameter vector
    μ = β₀ .+ β₁ * log.(flow) # find mean
    ll = sum(logpdf.(Normal.(μ, σ), tds)) # compute log-likelihood
    return ll
end

lb = [0.0, -1000.0, 1.0]
ub = [1000.0, 1000.0, 100.0]
θ₀ = [500.0, 0.0, 50.0]
optim_out = Optim.optimize(θ -> -tds_riverflow_loglik(θ, tds.tds_mgL, tds.discharge_cms), lb, ub, θ₀)
θ_mle = optim_out.minimizer; digits=0

x = 1:0.1:60
μ = θ_mle[1] .+ θ_mle[2] * log.(x)
p1 = scatter(tds.discharge_cms, tds.tds_mgL, markersize=3, xlabel=L"Log-Discharge (log-m$^3$/s)", ylabel="Total dissolved solids (mg/L)", label="Data", size=(500, 450))
plot!(p1, x, μ, linewidth=3, color=:red, label="Regression")

pred = θ_mle[1] .+ θ_mle[2] * log.(tds.discharge_cms)
resids = pred - tds.tds_mgL
p2 = scatter(log.(tds.discharge_cms), resids, markersize=3, ylabel=L"$r$ (mg/L)", xlabel=L"Log-Discharge (log-m$^3$/s) ", size=(500, 450))
p3 = scatter(resids[1:end-1], resids[2:end], markersize=3, ylabel=L"$r_{t+1}$ (mg/L)", xlabel=L"$r_t$ (mg/L)", size=(500, 450))
T = nrow(tds)
rlr = lm([ones(T-1) resids[1:end-1]], resids[2:end])
rpred = predict(rlr, [ones(T-1) resids[1:end-1]])
plot!(p3, resids[1:end-1], rpred, linewidth=3, color=:red)

display(p1)
display(p2)
display(p3)
(a) Autocorrelation of linear regression residuals
(b)
(c)
Figure 7

Residuals have autocorrelation 0.13.

Autoregressive (AR) Models

AR(1) Models

AR(p): (autoregressive of order \(p\)):

\[ \begin{align*} y_t &= \alpha + \rho y_{t-1} + \varepsilon_t \\ \varepsilon &\sim \mathcal{D}(\theta) \\ y_0 &\sim \mathcal{G}(\psi) \end{align*} \]

\(\mathbb{E}[\varepsilon] = 0\), \(\text{cor}(\varepsilon_i, \varepsilon_j) = 0\), \(\text{cor}(\varepsilon_i, y_0) = 0\)

Uses of AR Models

AR models are commonly used for prediction: bond yields, prices, electricity demand, short-run weather.

But may have little explanatory power: what causes the autocorrelation?

AR(1) Models

(a) AR(1) model
(b)
Figure 8

Simulating AR(1) Data

  • Draw an initial \(y_0 \sim \mathcal{G}(\psi)\)
  • For \(t \geq 1\):
    • Draw \(\varepsilon_t \sim \mathcal{D}(\theta)\)
    • Set \(y_t = \alpha + \rho y_{t-1} + \varepsilon_t\)

Diagnosing Autocorrelation

Plot \(\varsigma(i)\) over a series of lags.

Data generated by an AR(1) with \(\rho = 0.7\).

Notice: No explicit higher order autocorrelation, but \(\varsigma(2) > 0\).

Figure 9: Autocorrelation Function

“Memory” of Autocorrelation

\[\begin{aligned} y_t &= \alpha + \rho y_{t-1} + \varepsilon_t \\ &= \alpha + \rho \left(\alpha + \rho y_{t-2} + \varepsilon_{t-1}\right) + \varepsilon_t \\ &= \alpha + \rho \alpha + \rho^2 y_{t-2} + \rho \varepsilon_{t-1} + \varepsilon_t \\ &= \ldots \\ &= \alpha \sum_{k=0}^{t-1} \rho^k + \rho^t y_0 + \sum_{k=0}^{t-1} \varepsilon_{t-k} \rho^k \end{aligned}\]

Partial Autocorrelation

Can isolate \(\varsigma(i)\) independent of \(\varsigma(i-k)\) through partial autocorrelation.

Typically estimated through regression, \[y = \sum_{k=1}^i \varsigma(i-k) y_{i-k}.\]

Figure 10: Partial autocorrelation Function

Key Points

Key Points

  • Time series exhibit serial dependence (autocorrelation)
  • AR(1): autoregressive model with explicit lag-1 dependence.
  • AR(1) is stationary if and only if \(|\rho| < 1\).
  • In general, AR models useful for forecasting, pretty useless for explanation.
  • Similar concepts in spatial data: spatial correlation (distance/adjacency vs. time), lots of different models.

For More On Time Series

Courses

  • STSCI 4550/5550 (Applied Time Series Analysis),
  • CEE 6790 (heavy focus on spectral analysis and signal processing)

Books

  • Shumway & Stoffer (2025)
  • Hyndman & Athanasopoulos (2021)
  • Durbin & Koopman (2012)
  • Banerjee et al. (2011)
  • Cressie & Wikle (2011)

Upcoming Schedule

Next Classes

Friday: AR(1) Inference

Next Week: Extreme Values (Theory and Models)

Next Unit: Hypothesis Testing, Model Evaluation, and Comparison

Assessments

HW2: Due Friday at 9pm.

Exercises: Available after class (will include GLMs + time series).

HW3: Assigned on 3/2, due 3/13.

References

References (Scroll For Full List)

Banerjee, S., Carlin, B. P., & Gelfand, A. E. (2011). Hierarchical modeling and analysis for spatial data, second edition (2nd ed.). Philadelphia, PA: Chapman & Hall/CRC. https://doi.org/10.1201/b17115
Cressie, N., & Wikle, C. K. (2011). Statistics for Spatio-Temporal Data. Hoboken, NJ: Wiley.
Durbin, J., & Koopman, S. J. (2012). Time series analysis by state space methods (2nd ed.). London, England: Oxford University Press. https://doi.org/10.1093/acprof:oso/9780199641178.001.0001
Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed). Melbourne, Australia: OTexts. Retrieved from https://otexts.com/fpp3/
Shumway, R. H., & Stoffer, D. S. (2025). Time series analysis and its applications. Springer.