Time series analysis is essential for analyzing data ordered chronologically,whether tracking stock prices, weather patterns, sales trends, or any sequential measurements. Understanding how to decompose time series, identify patterns, detect autocorrelation, and build predictive models enables you to extract actionable insights and make informed forecasts.

This comprehensive guide covers all aspects of time series analysis with practical R implementations.

Time Series Objects and Basic Operations

Creating Time Series Objects

# Create time series object
data <- c(100, 105, 110, 108, 115, 120, 118, 125)

# Create ts object (frequency = 1 means annual data)
ts_annual <- ts(data, start = 2015, frequency = 1)
print(ts_annual)

# Monthly data (frequency = 12)
monthly_data <- c(150, 152, 155, 158, 160, 162, 165, 168, 170, 172, 175, 178)
ts_monthly <- ts(monthly_data, start = c(2020, 1), frequency = 12)
print(ts_monthly)

# Quarterly data (frequency = 4)
quarterly_data <- c(100, 110, 105, 115, 120, 130, 125, 135)
ts_quarterly <- ts(quarterly_data, start = c(2020, 1), frequency = 4)
print(ts_quarterly)

# Access time series components
time(ts_monthly)       # Time index
frequency(ts_monthly)  # Frequency
start(ts_monthly)      # Start time
end(ts_monthly)        # End time

Visualizing Time Series

# Plot time series
plot(ts_monthly, main = "Monthly Data", xlab = "Time", ylab = "Value")

# Multiple series plot
plot(ts.union(ts_monthly, ts_monthly * 1.1), plot.type = "single",
     col = c("blue", "red"), main = "Comparison")

# Add grid and legend
plot(ts_monthly)
grid(TRUE)

Autocorrelation and Partial Autocorrelation

ACF (Autocorrelation Function)

ACF measures correlation between observations at different lags;essential for identifying patterns and dependencies.

# Calculate and plot autocorrelation
data <- c(12.39, 11.25, 12.15, 13.48, 13.78, 12.89, 12.21, 12.58)

# ACF plot
acf(data, main = "Autocorrelation Function")

# Get ACF values
acf_result <- acf(data, plot = FALSE)
print(acf_result)
# Lag 0: 1.000 (always 1)
# Lag 1: 0.396
# Lag 2: -0.407
# etc.

# Interpretation of ACF
# Blue dashed lines = confidence bounds
# Spikes beyond bounds = significant autocorrelation
# Gradual decay = non-stationary data
# Sharp cutoff = MA process

PACF (Partial Autocorrelation Function)

PACF shows correlation at lag k after removing effects of intermediate lags.

# Plot PACF
pacf(data, main = "Partial Autocorrelation Function")

# Get PACF values
pacf_result <- pacf(data, plot = FALSE)

# ACF vs PACF patterns:
# AR(p): PACF cuts off at lag p, ACF decays gradually
# MA(q): ACF cuts off at lag q, PACF decays gradually
# ARMA(p,q): Both decay gradually

Time Series Decomposition

Decomposition separates time series into components: trend, seasonality, and residuals.

# Built-in time series with trend and seasonality
data(AirPassengers)
head(AirPassengers)

# Additive decomposition (Trend + Seasonal + Residual)
decomposed_add <- decompose(AirPassengers, type = "additive")
plot(decomposed_add)

# Multiplicative decomposition (Trend × Seasonal × Residual)
decomposed_mult <- decompose(AirPassengers, type = "multiplicative")
plot(decomposed_mult)

# Extract components
trend_component <- decomposed_add$trend
seasonal_component <- decomposed_add$seasonal
residual_component <- decomposed_add$random

# Advanced: Using stl() for seasonal decomposition
stl_result <- stl(AirPassengers, s.window = "periodic")
plot(stl_result)

Interpreting Decomposition

# Trend: Long-term direction
plot(decomposed_add$trend, main = "Trend Component")

# Seasonality: Regular repeating patterns
plot(decomposed_add$seasonal, main = "Seasonal Component")

# Residuals: What remains after trend and seasonality
plot(decomposed_add$random, main = "Residual Component")

# Check for stationarity of residuals
acf(na.omit(decomposed_add$random))

Stationarity and Differencing

Stationary series have constant mean and variance;required for most time series models.

# Test for stationarity (Augmented Dickey-Fuller test)
library(tseries)

# Non-stationary series (random walk)
nonstationary <- cumsum(rnorm(100))
plot(nonstationary, main = "Non-stationary Series")

# ACF of non-stationary data (slow decay)
acf(nonstationary)

# Differencing to achieve stationarity
differenced <- diff(nonstationary)
plot(differenced, main = "Differenced Series")

# ACF of differenced data (quicker decay)
acf(differenced)

# Adf test
adf_test <- adf.test(differenced)
print(adf_test)  # p-value < 0.05 indicates stationarity

ARIMA Models

ARIMA(p,d,q): AutoRegressive Integrated Moving Average for forecasting stationary series.

# Fit ARIMA model
library(forecast)

data <- AirPassengers

# Auto selection of ARIMA parameters
model_auto <- auto.arima(data)
print(model_auto)
# Output: ARIMA(0,1,1)(0,1,1)[12]

# Manual ARIMA specification
# ARIMA(p,d,q)
# p: AR (autoregressive) order
# d: differencing (integration) order
# q: MA (moving average) order

model <- arima(data, order = c(0, 1, 1))
summary(model)

# Model coefficients and their significance
coefficients(model)

Forecasting

# Forecast using ARIMA model
library(forecast)

# Fit and forecast
model <- auto.arima(AirPassengers)
forecast_result <- forecast(model, h = 24)  # 24 months ahead

# Plot forecast
plot(forecast_result)

# Extract forecast values
point_forecast <- forecast_result$mean
upper_bound <- forecast_result$upper[, 2]  # 95% upper
lower_bound <- forecast_result$lower[, 2]  # 95% lower

# Create forecast data frame
forecast_df <- data.frame(
  Month = seq(length(AirPassengers) + 1, length(AirPassengers) + 24),
  Forecast = as.numeric(point_forecast),
  Lower = as.numeric(lower_bound),
  Upper = as.numeric(upper_bound)
)

# Exponential smoothing forecast
model_ets <- ets(AirPassengers)
forecast_ets <- forecast(model_ets, h = 24)
plot(forecast_ets)

Seasonality Detection and Handling

# Seasonal decomposition
seasonality <- decomposed_add$seasonal
plot(seasonality, main = "Seasonal Pattern")

# Seasonal adjustment (remove seasonality)
seasonally_adjusted <- AirPassengers - decomposed_add$seasonal
plot(seasonally_adjusted, main = "Seasonally Adjusted Data")

# SARIMA for seasonal ARIMA
# SARIMA(p,d,q)(P,D,Q)s includes seasonal components
model_sarima <- arima(AirPassengers, order = c(1, 1, 1),
                      seasonal = list(order = c(1, 1, 1), period = 12))
summary(model_sarima)

# Forecast with SARIMA
forecast_sarima <- forecast(model_sarima, h = 24)
plot(forecast_sarima)

Model Diagnostics

# Residual analysis
model <- auto.arima(AirPassengers)

# Residuals should be:
# 1. White noise (no autocorrelation)
# 2. Normally distributed
# 3. Mean zero

residuals <- residuals(model)

# Check residuals
par(mfrow = c(2, 2))
plot(residuals, main = "Residuals Over Time")
acf(residuals, main = "ACF of Residuals")
pacf(residuals, main = "PACF of Residuals")
hist(residuals, main = "Distribution of Residuals")
par(mfrow = c(1, 1))

# Formal tests
# Ljung-Box test (autocorrelation)
Box.test(residuals, lag = 10, type = "Ljung-Box")

# Shapiro-Wilk test (normality)
shapiro.test(residuals)

Practical Time Series Workflow

# Complete workflow example
print("=== Time Series Analysis Workflow ===")

# 1. Load and visualize data
data(AirPassengers)
plot(AirPassengers)

# 2. Check stationarity
library(tseries)
adf <- adf.test(AirPassengers)
cat("ADF test p-value:", adf$p.value, "\n")

# 3. Decompose
decomposed <- decompose(AirPassengers)
plot(decomposed)

# 4. Difference if needed
differenced <- diff(AirPassengers)

# 5. Fit model
library(forecast)
model <- auto.arima(AirPassengers)
summary(model)

# 6. Check diagnostics
checkresiduals(model)

# 7. Forecast
forecast_result <- forecast(model, h = 24)
plot(forecast_result)

# 8. Accuracy
accuracy(forecast_result)

Best Practices

  1. Visualize first - Always plot your time series before analysis
  2. Check stationarity - Use ADF test or visual inspection
  3. Remove seasonality - Use differencing or seasonal decomposition
  4. Check ACF/PACF - Guide ARIMA parameter selection
  5. Validate model - Use holdout test set for forecast evaluation
  6. Monitor residuals - Should be white noise
  7. Consider external factors - Events, holidays, structural breaks
  8. Use ensemble methods - Combine multiple forecasts for better accuracy

Common Questions

Q: What does ACF cutoff at lag 2 mean? A: MA(2) process - moving average depends on past 2 errors

Q: Should I difference monthly data with seasonality? A: Yes, use SARIMA with seasonal differencing (D=1)

Q: How do I choose between ARIMA and exponential smoothing? A: Try both with auto.arima() and ets(), compare on test set

Q: What if my data is not stationary? A: Difference it (d=1 or d=2 in ARIMA) until stationary

Q: How many forecasts ahead can I make? A: Generally 1-2 seasons ahead is reliable; longer forecast = wider confidence intervals

Download R Script

Get all code examples from this tutorial: time-series-examples.R