SMHS TimeSeries

Scientific Methods for Health Sciences - Time Series Analysis

• See the DSPA2 Extensible Longitudinal Data Analysis Chapter.
• The SOCR Statistical Data Analyzer (SDA) provides live app-based data and multiple tools for interactive time-varying data analytics.

Overview

Time series data is a sequence of data points measured at successive, equally spaced points in time. Time series analysis is utilized across various disciplines, such as monitoring industrial processes, tracking business metrics, and analyzing longitudinal health records. In this module, we present a comprehensive introduction to time series data and explore both classical and modern techniques used in this rapidly growing field. The goal is to extract meaningful statistics, identify underlying patterns (such as trends and seasonality), and construct predictive models. We will illustrate the application of these techniques using interactive examples in R.

Motivation

Economic data (e.g., daily share prices, monthly sales, annual income) and biophysical data (e.g., daily temperature, ECG readings, epidemiological counts) are all examples of time series data. A time series is an ordered sequence of values of a variable measured at equally spaced time intervals. What are the most effective ways to model time series data? How can we extract actionable information, make inferences about underlying mechanisms, and forecast future observations? Answering these questions requires understanding the mathematical properties of time series processes and applying robust statistical modeling techniques.

Theory

Components of Time Series

A time series \(x_t\) can generally be decomposed into three primary components, plus an irregular error term:

• Trend Component (\(T_t\)): A long-run increase or decrease over time (overall upward or downward movement). The trend can be linear or non-linear.
• Seasonal Component (\(S_t\)): Short-term, regular wave-like patterns completing a cycle within a fixed period (e.g., daily, monthly, quarterly).
• Cyclical Component (\(C_t\)): Long-term wave-like patterns with variable lengths, often measured peak-to-peak or trough-to-trough, typically tied to economic or biological cycles.
• Irregular Component (\(I_t\)): Unpredictable, random residuals, or "noise," which may be due to random variations or unusual events.

These components can combine additively or multiplicatively:

• Additive Model: \(x_{t}=T_{t}+S_{t}+C_{t}+I_{t}\)
• Multiplicative Model: \(x_{t}=T_{t} \times S_{t} \times C_{t} \times I_{t}\)

A logarithmic transformation can convert a multiplicative model into an additive one\[\log(x_{t})=\log(T_{t})+\log(S_{t})+\log(C_{t})+\log(I_{t})\].

• Probabilistic models: The components still add or multiply together across time.

Stationarity

Most statistical forecasting methods assume that a time series is approximately stationary. A time series is said to be strictly stationary if its statistical properties are unaffected by a shift in time. In practice, we rely on weak stationarity, which requires: 1. The mean is constant over time\[E[x_t] = \mu\] for all \(t\). 2. The variance is finite and constant\[Var(x_t) = \gamma(0) < \infty\]. 3. The autocovariance between any two points depends only on the time lag \(h\) between them\[Cov(x_t, x_{t+h}) = \gamma(h)\] for all \(t\) and \(h\).

Real-world time series are usually non-stationary (exhibiting trends or varying variance) but can often be transformed to achieve stationarity (e.g., differencing, log transformations).

Autocorrelation and Partial Autocorrelation

To characterize a time series, we examine its mean and covariance structure. For a stationary process, the Autocovariance Function (ACVF) at lag \(h\) is defined as\[\gamma(h) = E[(x_t - \mu)(x_{t+h} - \mu)]\] The Autocorrelation Function (ACF) standardizes this by the variance \(\gamma(0)\)\[\rho(h) = \frac{\gamma(h)}{\gamma(0)}\] Properties of the ACF include\[\rho(0) = 1\], \(|\rho(h)| \le 1\), and \(\rho(h) = \rho(-h)\) (it is an even function).

The Partial Autocorrelation Function (PACF) measures the correlation between \(x_t\) and \(x_{t+h}\) after removing the linear dependence on the intermediate observations \(x_{t+1}, \dots, x_{t+h-1}\). The ACF and PACF are crucial tools for identifying the order of Autoregressive (AR) and Moving Average (MA) models.

Estimation of ACF

Given \(T\) observations, the sample autocovariance is estimated as\[\hat{\gamma}(h) = \frac{1}{T} \sum_{t=1}^{T-h} (x_t - \bar{x})(x_{t+h} - \bar{x})\] The sample ACF is \(\hat{\rho}(h) = \hat{\gamma}(h) / \hat{\gamma}(0)\). Under the null hypothesis that the true process is White Noise, the large-sample standard error of the sample ACF is approximately\[SE(\hat{\rho}(h)) \approx \frac{1}{\sqrt{T}}\] A common rule of thumb is that 95% of the sample autocorrelations for a White Noise process should fall within the bounds \(\pm 2/\sqrt{T}\).

White Noise

A time series \(w_t\) is White Noise (\(WN\)) if it is weakly stationary with\[\rho(h) = \begin{cases} 1 & \text{if } h=0 \\ 0 & \text{if } h \neq 0 \end{cases}\] A special case is Gaussian White Noise, denoted \(w_t \sim \text{WN}(0, \sigma^2)\), where \(w_t\) is independently and identically distributed (i.i.d.) according to a normal distribution.

Backshift Notation and Linear Operators

Modern time series theory heavily utilizes the backshift operator \(B\), defined as \(Bx_t = x_{t-1}\).

• Differencing can be written as \(\nabla x_t = (1-B)x_t\).
• An AR(1) model \(x_t = \phi x_{t-1} + w_t\) is written as \((1-\phi B)x_t = w_t\).

Autoregressive (AR) Models

An autoregressive model of order \(p\), denoted AR(\(p\)), models the current value as a linear combination of its past values plus a shock\[x_t = \phi_1 x_{t-1} + \phi_2 x_{t-2} + \dots + \phi_p x_{t-p} + w_t\] Using backshift notation\[\phi(B)x_t = w_t\], where \(\phi(B) = 1 - \phi_1 B - \dots - \phi_p B^p\). Assumptions: The process must be stationary, which requires that the roots of \(\phi(B) = 0\) lie strictly outside the unit circle. ACF/PACF Patterns: The ACF of a stationary AR(\(p\)) process decays gradually (tails off), while the PACF cuts off abruptly after lag \(p\).

Moving Average (MA) Models

A moving average model of order \(q\), denoted MA(\(q\)), models the current value as a linear combination of current and past white noise shocks\[x_t = w_t + \theta_1 w_{t-1} + \theta_2 w_{t-2} + \dots + \theta_q w_{t-q}\] Using backshift notation\[x_t = \theta(B)w_t\], where \(\theta(B) = 1 + \theta_1 B + \dots + \theta_q B^q\). Assumptions: MA processes are always stationary. However, they must be invertible, meaning the roots of \(\theta(B) = 0\) lie strictly outside the unit circle. Invertibility ensures the process can be represented as an infinite-order AR model. ACF/PACF Patterns: The ACF of an MA(\(q\)) process cuts off abruptly after lag \(q\), while the PACF decays gradually (tails off).

ARIMA Models

The Autoregressive Integrated Moving Average model, ARIMA(\(p, d, q\)), combines AR and MA models after differencing the data \(d\) times to achieve stationarity\[\phi(B)(1-B)^d x_t = \theta(B)w_t\] If seasonality is present, we extend this to the Seasonal ARIMA (SARIMA) model, which incorporates seasonal AR, MA, and differencing terms.

In practice, R's forecast package fits and reports ARIMA models with seasonality, e.g.,

ARIMA(p,d,q)(P,D,Q)[m]

These six parameters define how the model handles trends, cycles, and shocks in the data.

For instance, this model ARIMA(2,1,1)(0,1,0)[12] has two parts, the Non-Seasonal component (first 3 parameters), and the Seasonal component, (the last 3 parameters) .

• Non-Seasonal Component: (p, d, q): These parameters handle the short-term relationships between consecutive observations.
  • p (Autoregressive - AR): The number of lag observations included in the model. A value of 2 means the model uses the two previous time points to predict the next one.
  • d (Degree of Differencing): The number of times the raw observations are differenced to make the data stationary (removing trends). A value of 1 means the model is looking at the change between points rather than the absolute values.
  • q (Moving Average - MA): The size of the moving average window applied to forecast errors. A value of 1 means the model accounts for the error made in the previous prediction.

• Seasonal Component: (P, D, Q): These parameters handle repeating patterns (like the AirPassengers monthly spikes).
  • P (Seasonal Autoregressive): Similar to p, but specifically for the seasonal lags. For monthly data, P=1 would mean using the value from the same month last year to predict this month, in the output, this is 0.
  • D (Seasonal Differencing): The number of seasonal differences. A value of 1 means the model subtracted last year's value from this year's value to remove seasonal trends.
  • Q (Seasonal Moving Average): Similar to q, but for seasonal forecast errors, in the output, this is 0.
  • [m] (Seasonal Period): The number of periods in each season. For the AirPassengers dataset, this is 12 because the data is monthly and repeats every year.

For a specific model: ARIMA(2,1,1)(0,1,0)[12]

Parameter	Value	Meaning
p	2	Uses 2 lags for the AR part.
d	1	Applied 1st order differencing to remove the general trend.
q	1	Uses 1 lag for the Moving Average error correction.
P	0	No seasonal autoregressive terms.
D	1	Applied seasonal differencing (Value - Value same month last year).
Q	0	No seasonal moving average terms.
m	12	Data follows a 12-month annual cycle.

• Control Parameters: The Stepwise and Approximation Flags, stepwise=FALSE and approximation=FALSE, force auto.arima to search more thoroughly through all possible combinations of these 6 parameters rather than using a shortcut search. This results in a better fit but takes longer to compute.

Methods and Data Analysis

Exploratory Data Analysis (EDA) and Diagnostics

Before modeling, we must visualize the data, check for outliers, and assess distributional assumptions. A standard approach involves time series plots, Q-Q plots (to check normality), and examining residuals.

If we fit a simple mean model \(x_t = \mu + w_t\), the residuals \(e_t = x_t - \hat{\mu}\) should resemble White Noise. We diagnose this by checking: 1. Zero mean of residuals. 2. Constant variance across time (homoscedasticity). 3. No autocorrelation (checked via ACF plots of residuals). 4. Normality of residuals (checked via Q-Q plots).

Smoothing Techniques

Smoothing helps identify underlying structures by averaging out the noise.

Moving Averages

A simple centered moving average smooths \(x_t\) using a window of size \(k\)\[v_t = \frac{1}{k} \sum_{j=-(k-1)/2}^{(k-1)/2} x_{t+j}\]

# Simulate White Noise and apply Moving Averages
set.seed(123)
w <- ts(rnorm(150))  # generate 150 data points from standard normal distribution
ma3 <- filter(w, sides=2, rep(1/3, 3))  # 3-period centered moving average
ma9 <- filter(w, sides=2, rep(1/9, 9))  # 9-period centered moving average

plot.ts(w, main="White Noise vs. Moving Averages", ylab="Value")
lines(ma3, col="red", lwd=2)
lines(ma9, col="blue", lwd=2)
legend("topright", legend=c("WN", "MA(3)", "MA(9)"), col=c("black", "red", "blue"), lwd=2)

• Note: A larger window (lag) results in a smoother time series plot, but introduces more lag in detecting changes.*

Exponential Smoothing

Unlike moving averages, exponential smoothing assigns exponentially decreasing weights as observations get older.

• Simple Exponential Smoothing (SES): Appropriate for data with no trend or seasonality. The forecast is \(\hat{x}_{t+1} = \alpha x_t + (1-\alpha)\hat{x}_t\), where \(0 < \alpha < 1\) is the smoothing parameter.
• Holt-Winters Method: Extends SES to capture trends (additive or multiplicative) and seasonality.

# Using built-in AirPassengers dataset
data(AirPassengers)
plot(AirPassengers, main="AirPassengers Dataset", ylab="Passengers")

# Holt-Winters Additive Model
hw_add <- HoltWinters(AirPassengers, seasonal="additive")
plot(hw_add, main="Holt-Winters Additive Filtering")

Nonparametric Smoothing

• Loess (Local Polynomial Regression): Fits simple models to localized subsets of data. Controlled by the `span` parameter. A larger span yields a smoother curve.
• Splines: Piecewise polynomials joined smoothly at "knots." Controlled by degrees of freedom or a smoothing parameter (`spar`).

# Loess and Spline Smoothing on a time series
t <- time(AirPassengers)
fit_loess <- lowess(t, AirPassengers, f=0.3) # f is the span
fit_spline <- smooth.spline(t, AirPassengers, spar=0.8)

plot(AirPassengers, main="Nonparametric Smoothing", ylab="Passengers")
lines(fit_loess, col="red", lwd=2)
lines(fit_spline, col="blue", lwd=2, lty=2)
legend("topleft", legend=c("Loess", "Spline"), col=c("red", "blue"), lwd=2, lty=c(1,2))

Modeling Seasonality and Trend

Classical Decomposition

We can decompose a time series into its trend, seasonal, and irregular components using moving averages or Local Regression (LOESS).

# Multiplicative decomposition using LOESS
decomp <- stl(AirPassengers, s.window="periodic")
plot(decomp, main="Classical Decomposition via STL")

Regression with Dummy Variables

Seasonality can be explicitly modeled using regression with dummy variables. For monthly data with an intercept, we use 11 dummy variables to avoid the dummy variable trap.

# Regression with seasonal dummy variables
y <- log(AirPassengers) # Log transform to stabilize variance
t <- time(y)
Q <- factor(cycle(y))   # Creates monthly factors (1-12)

# Fit linear regression with trend and seasonality
reg_model <- lm(y ~ t + Q)
summary(reg_model)

# Diagnostics of the regression model
par(mfrow=c(2,2))
plot(reg_model, main="Regression Diagnostics")

Modern ARIMA Modeling (Box-Jenkins Approach)

The modern standard for time series forecasting involves the Box-Jenkins iterative cycle: Identification, Estimation, and Diagnostic Checking.

1. Identification: Use ACF and PACF plots, and unit root tests (like Augmented Dickey-Fuller) to determine \(d\). Use algorithms to suggest \(p\) and \(q\).

2. Estimation: Maximize the likelihood function to estimate AR and MA coefficients.

3. Diagnostic Checking: Ensure residuals are White Noise.

The forecast package in R automates this efficiently:

# Install and load forecast package if necessary
# install.packages("forecast")
library(forecast)

# dev.off()  # reset display, if necessary

# auto.arima automatically identifies the best SARIMA model
fit <- auto.arima(AirPassengers, seasonal=TRUE, stepwise=FALSE, approximation=FALSE)
summary(fit)

# Check residuals to ensure they are White Noise
checkresiduals(fit)

# Forecast the next 24 months
fc <- forecast(fit, h=24)
plot(fc, main="24-Month Forecast using SARIMA")

Next, let's try to fit manually another (specific) ARIMA(1,1,5) model where the (p,d,q) parameters are hard-coded. Then, compare fitManual to the optimized model (fit <- auto.arima()), by contrasting the corresponding AIC, BIC, and log-likelihood estimates of each model, smaller values indicate better model and higher model fidelity.

Also, compare the 24-month forward forecasts of fitManual vs. the optimized model fit.

fitManual <- arima(AirPassengers, order=c(1,1,5)); summary(fitManual)

fcManual <- forecast(fitManual, h=24)
plot(fcManual, main="(Manual) 24-Month Forecast using SARIMA")

fcManual$model$aic  # [1] 1353.116
fcManual$model$loglik  # [1] -669.5579

Outlier Detection

Outliers can severely distort time series models (inflating variance, skewing ACF). They can be detected visually via time plots and boxplots, or systematically using residuals.

# Creating a dataset with injected outliers
set.seed(42)
dat <- ts(c(28,22,36,26,28,28,26,24,32,30,27,24,33,21,36,32,31,25,24,25,28,36,27,32,34,30,
            25,26,26,25,-44,23,21,30,33,29,27,29,28,22,26,27,16,31,29,36,32,28,40,19,
            37,23,32,29,-2,24,25,27,24,16,29,20,28,27,39,23))

par(mfrow=c(1,2))
plot.ts(dat, main="Time Series with Outliers")
boxplot(dat, main="Boxplot Identifying Outliers")

# Identifying and removing outliers based on extreme boxplot bounds
outlier_bounds <- boxplot.stats(dat)$out
dat_clean <- dat[!(dat %in% outlier_bounds)]

# Re-plot clean data
par(mfrow=c(1,1))
plot.ts(dat_clean, main="Cleaned Time Series", col="blue", lwd=2)

Applications

1. Epidemiology and Public Health: Time series models are essential for tracking disease outbreaks (e.g., COVID-19, influenza). SARIMA models help forecast case counts, accounting for weekly seasonal cycles (e.g., higher reporting on weekdays) and long-term trends, thereby informing hospital resource allocation. 2. Health Economics: Forecasting hospital readmission rates, healthcare costs, and resource utilization over time using ARIMA and Exponential Smoothing to optimize budget planning. 3. Clinical Monitoring: Analyzing ECG readings or continuous glucose monitoring data. Signal processing techniques combined with MA filtering help isolate true physiological signals from high-frequency noise. 4. Macroeconomics: As investigated by Nelson and Plosser (1982), understanding whether macroeconomic time series (like GDP or stock prices) are stationary around a deterministic trend or non-stationary stochastic processes fundamentally changes how we interpret economic shocks and formulate policy.

Software

Modern time series analysis in R relies on a few key packages:

• stats (Base R): Contains foundational functions like `ts()`, `acf()`, `pacf()`, `arima()`, `HoltWinters()`, and `stl()`.
• forecast (now transitioning to fable): Provides `auto.arima()`, `checkresiduals()`, `forecast()`, and modern tidy time series frameworks (`tsibble`).
• astsa: Applied Statistical Time Series Analysis, highly recommended for its companion datasets and straightforward functions like `acf2()` (combines ACF and PACF in one plot).
• tseries: Provides unit root tests like `adf.test()` to check for stationarity.

Problems

1) Consider a signal-plus-noise model of the general form \(x_{t}=s_{t}+w_{t}\), where \(w_{t}\) is Gaussian white noise with \(\sigma^{2}_{w} = 1\). Simulate and plot \(n=200\) observations from each of the following two models.

(a) \(x_{t} = s_{t}+w_{t}\), for \(t=1,2,\dots,200\), where \(s_{t}=\begin{cases} 0 & \text{ if } t=1,2,\dots,100 \\ 10\exp \left\{ -(t-100)/20 \right\} \cos\left( 2\pi t/4 \right) & \text{ if } t=101,102,\dots,200 \end{cases}\)

(b) \(x_{t}=s_{t}+w_{t}\), for \(t=1,2,\dots,200\), where \(s_{t}=\begin{cases} 0 & \text{ if } t=1,2,\dots,100 \\ 10\exp \left\{ -(t-100)/200 \right\} \cos\left( 2\pi t/4 \right) & \text{ if } t=101,102,\dots,200 \end{cases}\)

(c) Compare the signal modulators \((a) \exp \left\{-t /20 \right\}\) and \((b) \exp \left\{-t/200 \right\}\), for \(t=1,2,\dots,100\). What is the primary difference in the persistence of the signal?

2) (a) Generate \(n=100\) observations from the autoregression \(x_{t}=-0.9 x_{t-2} + w_{t}\), with \(\sigma_w=1\). Next, apply the moving average filter \(v_{t}=\left(x_{t}+x_{t-1}+x_{t-2}+x_{t-3} \right)/4\) to the data you generated. Now, plot \(x_{t}\) as a line and superimpose \(v_{t}\) as a dashed line. Comment on the behavior of \(x_{t}\) and how applying the moving average filter changes that behavior. [Hint: use v <- filter(x, rep(1/4,4), sides=1) for the filter].

(b) Repeat (a) but with \(x_{t}=\cos \left( 2 \pi t/4 \right)\).

(c) Repeat (b) but with added \(N(0,1)\) noise, \(x_{t}=\cos \left(2 \pi t/4 \right) + w_{t}\).

(d) Compare and contrast the smoothing effects observed in \((a) - (c)\).

3) For the two series, \(x_{t}\) in Problem 1 (a) and (b):

(a) Compute and plot the mean functions \(\mu_{x}(t)\) for \(t=1,2,\dots,200\).

(b) Explain why these models are not stationary.

4) Consider the time series \(x_{t} = \beta_{1} + \beta_{2} t + w_{t}\), where \(\beta_{1}\) and \(\beta_{2}\) are known constants and \(w_{t}\) is a white noise process with variance \(\sigma^{2}_{w}\).

(a) Determine whether \(x_{t}\) is stationary. Justify your answer mathematically.

(b) Show that the process \(y_{t}=x_{t} - x_{t-1}\) is stationary.

(c) Show that the mean of the moving average \(v_{t}= \frac{1}{2q+1} \sum_{j=-q}^{q} x_{t-j}\) is \(\beta_{1} + \beta_{2} t\), and provide a simplified expression for its autocovariance function.

5) A time series with a periodic component can be constructed from \(x_{t} = U_{1} \sin(2 \pi \omega_{0} t) + U_{2} \cos(2 \pi \omega_{0} t)\), where \(U_{1}\) and \(U_{2}\) are independent random variables with zero means and \(E(U^{2}_{1}) = E(U^{2}_{2}) = \sigma^{2}\). The constant \(\omega_{0}\) determines the period. Show that this series is weakly stationary with autocovariance function \(\gamma(h) = \sigma^{2} \cos(2 \pi \omega_{0} h)\). [Hint: you will need to refer to standard trigonometric identities for products of sine and cosine functions].

6) Suppose we would like to predict a single stationary series \(x_{t}\) with zero mean and autocorrelation function \(\gamma(h)\) at some time in the future, say \(t+l\), for \(l>0\).

(a) If we predict using only \(x_{t}\) and some scale multiple \(A\), show that the mean-square prediction error \(MSE(A)=E \left[ (x_{t+l}-A x_{t})^{2} \right]\) is minimized by the value \(A = \rho(l)\).

(b) Show that the minimum mean-square prediction error is \(MSE(A)= \gamma(0) \left[ 1- \rho^{2}(l) \right]\).

(c) Show that if \(x_{t+l} = Ax_{t}\) perfectly, then \(\rho(l) = 1\) if \(A > 0\), and \(\rho(l) = -1\) if \(A < 0\).

7) Let \(w_{t}\), for \(t=0, \pm 1, \pm 2, \dots\) be a normal white noise process, and consider the series \(x_{t}=w_{t} w_{t-1}\). Determine the mean and autocovariance function of \(x_{t}\), and state whether it is strictly stationary, weakly stationary, or neither.

8) Suppose \(x_{1}, x_{2}, \dots, x_{n}\) is a sample from the process \(x_{t} = \mu + w_{t} - 0.8 w_{t-1}\), where \(w_{t} \sim \text{WN}(0, \sigma^{2}_{w})\).

(a) Show that the mean function is \(E(x_{t})= \mu\).

(b) Calculate the standard error of \(\bar{x}\) for estimating \(\mu\).

(c) Compare your result in (b) to the case where \(x_{t}\) is pure white noise. Is the standard error smaller or larger? Explain the result in the context of positive/negative autocorrelation.

9) Although the model in Problem 1 (a) is not stationary, the sample ACF can still be computed and informative. For the data you generated in that problem, calculate and plot the sample ACF. Why does the ACF fail to decay to zero rapidly?

10)

(a) Simulate a series of \(n=500\) Gaussian white noise observations and compute the sample ACF, \(\hat{\rho}(h)\), to lag 20. Compare the sample ACF you obtain to the theoretical ACF \(\rho(h)\). Do approximately 95% of the sample autocorrelations fall inside the bounds \(\pm 2/\sqrt{500}\)?

(b) Repeat part (a) using only \(n=50\). How does changing \(n\) affect the variability of the sample ACF and the reliability of the confidence bounds?

References

• Cryer, J. D., & Chan, K. S. (2008). Time Series Analysis: With Applications in R. Springer.
• Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed.). OTexts. (Available freely online at otexts.com/fpp3)
• Nelson, C. R., & Plosser, C. I. (1982). Trends and random walks in macroeconomic time series: Some evidence and implications. Journal of Monetary Economics, 10(2), 139-162.
• Shumway, R. H., & Stoffer, D. S. (2017). Time Series Analysis and Its Applications: With R Examples (4th ed.). Springer.

---

• SOCR Home page: https://socr.umich.edu

---

---

Translate this page:

(default)	Deutsch	Español	Français	Italiano	Português	日本語	България	الامارات العربية المتحدة	Suomi	इस भाषा में	Norge
한국어	中文	繁体中文	Русский	Nederlands	Ελληνικά	Hrvatska	Česká republika	Danmark	Polska	România	Sverige