pacman::p_load(forecast, imputeTS, trend, tidyverse, gridExtra, tseries)

theme_set(theme_bw(base_size = 12))

Introduction

This lesson demonstrates a comprehensive time series analysis of the built-in AirPassengers dataset, containing monthly totals of international airline passengers from 1949 to 1960. We’ll cover:

Data exploration and visualization
Missing value imputation techniques
Time series decomposition
Trend analysis with Sen’s slope
Stationarity testing
Model fitting and forecasting

Data Preparation

Load and Examine Data

data("AirPassengers")
ap = AirPassengers

# Basic information
class(ap)

## [1] "ts"

start(ap)

## [1] 1949    1

end(ap)

## [1] 1960   12

frequency(ap)

## [1] 12

Initial Visualization

autoplot(ap) + 
  ggtitle("Monthly Airline Passengers (1949-1960)") +
  ylab("Passengers (thousands)") +
  xlab("Year") +
  theme_minimal()

Interpretation:

- Clear upward trend over time

- Increasing seasonal amplitude (multiplicative seasonality)

- No obvious missing values in this complete dataset (we’ll simulate missing values for demonstration)

Missing Value Imputation

Simulate Missing Values (5% random)

set.seed(123)
ap_with_na = ap
na_indices = sample(length(ap), size = 0.05*length(ap))
ap_with_na[na_indices] = NA

# Visualize missing values
ggplot_na_distribution(ap_with_na) +
  ggtitle("Distribution of Simulated Missing Values")

Imputation Methods Comparison

# Mean imputation
ap_mean = na_mean(ap_with_na)

# Linear interpolation
ap_linear = na_interpolation(ap_with_na)

# Seasonal decomposition imputation
ap_seadec = na_seadec(ap_with_na, algorithm = "interpolation")

# Kalman smoothing
ap_kalman = na_kalman(ap_with_na)

# Create comparison plot
plots = list()
plots[[1]] = ggplot_na_imputations(ap_with_na, ap_mean) + ggtitle("a. Mean Imputation")
plots[[2]] = ggplot_na_imputations(ap_with_na, ap_linear) + ggtitle("b. Linear Interpolation")
plots[[3]] = ggplot_na_imputations(ap_with_na, ap_seadec) + ggtitle("c. Seasonal Decomposition")
plots[[4]] = ggplot_na_imputations(ap_with_na, ap_kalman) + ggtitle("d. Kalman Smoothing")

grid.arrange(grobs = plots, ncol = 2)

Interpretation:

- Mean imputation performs poorly for time series data

- Linear interpolation works better but ignores seasonality

- Seasonal decomposition and Kalman filtering produce the most realistic imputations

- We’ll proceed with the seasonally adjusted imputation

Time Series Decomposition

ap_stl = stl(ap_kalman, s.window = "periodic", robust = TRUE)

autoplot(ap_stl) +
  ggtitle("STL Decomposition with Robust Weights")

Components:

- Trend: Clear upward trajectory accelerating over time

- Seasonal: Consistent yearly pattern with peak in summer

- Random: Residuals show some pattern, suggesting limitations of classical decomposition

Advantages over stl:

- Handles changing seasonal patterns better

- More robust to outliers

- Allows for non-periodic seasonality

Trend Analysis

Sen’s Slope Estimator

# Convert to annual sums for trend analysis
ap_yearly = aggregate(ap_kalman, nfrequency = 1, FUN = sum)

# Perform Mann-Kendall test and Sen's slope
mk_test = mk.test(as.numeric(ap_yearly))
sens_slope = sens.slope(as.numeric(ap_yearly))

cat("Mann-Kendall Test p-value:", mk_test$p.value, "\n")

## Mann-Kendall Test p-value: 8.303107e-06

cat("Sen's Slope Estimate:", sens_slope$estimates, "passengers/year\n")

## Sen's Slope Estimate: 378.056 passengers/year

cat("Confidence Interval:", sens_slope$conf.int, "\n")

## Confidence Interval: 345.2036 413.3217

Interpretation:

- Highly significant upward trend (p < 0.001)

- Estimated annual increase of 378 thousand passengers

- 95% CI: 345 to 413 thousand passengers/year

Trend Visualization

ggplot(data.frame(Time = time(ap_yearly), Passengers = as.numeric(ap_yearly)),
       aes(x = Time, y = Passengers)) +
  geom_line(color = "steelblue", size = 1) +
  geom_smooth(method = "lm", se = FALSE, color = "red", linetype = "dashed") +
  labs(title = "Annual Passenger Totals with Linear Trend",
       y = "Passengers (thousands)")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
## `geom_smooth()` using formula = 'y ~ x'

Stationarity Analysis

Augmented Dickey-Fuller Test

adf_test = adf.test(ap_kalman)

## Warning in adf.test(ap_kalman): p-value smaller than printed p-value

cat("ADF p-value:", adf_test$p.value, "\n")

## ADF p-value: 0.01

# Test after differencing
adf_diff = adf.test(diff(ap_kalman))

## Warning in adf.test(diff(ap_kalman)): p-value smaller than printed p-value

cat("ADF after differencing p-value:", adf_diff$p.value, "\n")

## ADF after differencing p-value: 0.01

Interpretation:

Original series is stationary (p = 0.01), which is unusual for most time series
Differenced series is also stationary (p = 0.01)
Use ARMA (AutoRegressive Moving Average, d = 0) for stationary and ARIMA (AutoRegressive Moving Average) for non-stationary

Seasonal Differencing

Seasonal differencing means subtracting the values from the corresponding month of the previous year. This makes the series stationary, meaning it has constant mean and variance over time.

ggtsdisplay(diff(ap_kalman, lag = 12), 
            main = "Seasonally Differenced Series")

Observations:

ACF - Autocorrelation function, PACF - Partial autocorrelation function
Autocorrelation means the relationship between a data point and its past values at specific time intervals (lag)
Sesonal difference removes most of the seasonal pattern
Some remaining autocorrelation suggests need for further modeling

Model Building

SARIMA Model Identification (Seasonal AutoRegressive Integrated Moving Average)

ACF and PACF plots help identify the order of the AR (AutoRegressive) and MA (Moving Average) components in a SARIMA model. These require identifying the parameters p, d, q for the non-seasonal part and P, D, Q for the seasonal part. We need stationary series for ACF and PACF analysis, so we will use the seasonally differenced series. Lags: ACF or PACF correlations are shown based on different lags. The significant correlations are outside the 95% confidence interval (blue dashed lines).

ACF/PACF in original series: a gradual decay suggests an AR process, while a sharp cut-off suggests an MA process. ACF/PACF in original series: a sharp cut-off suggests an AR process, while a gradual decay suggests an MA process.

The spikes in the ACF and PACF plots indicate the presence of autocorrelation at specific lags, which can help determine the order of the AR and MA components in the SARIMA model. For example, if the ACF cuts off after lag 1, it suggests an MA(1) model, while if the PACF cuts off after lag 1, it suggests an AR(1) model.

ARMA model is applied when ACF and PACF both show gradual decay, indicating a combination of AR and MA processes.

Seasonality: Signifiant spikes at seasonal lags (e.g., 12, 24) suggest seasonal AR or MA terms.

For ARIMA models, we need to identify the order of differencing (d) required to make the series stationary. The ACF and PACF plots can help determine the appropriate differencing order.

Differencing order means the number of times we subtract the previous value from the current value to make the series stationary. For example, if we need to subtract the previous value once, we use d = 1.

p, q, d, P, Q, D are the parameters of the SARIMA model, where: p = order of the AR component, q = order of the MA component, d = degree of differencing, P = order of the seasonal AR component, Q = order of the seasonal MA component, D = degree of seasonal differencing.

Therefore, ACF and PACF plots are essential for identifying the appropriate parameters for the SARIMA model, which can then be used for forecasting future values in the time series.

# Original series
p1 = ggAcf(ap_kalman) + ggtitle("ACF of Original Series") 
p2 = ggPacf(ap_kalman) + ggtitle("PACF of Original Series") 

# Differenced series with suitably lagged and iterated differences
p3 = ggAcf(diff(ap_kalman)) + ggtitle("ACF of Differenced Series") 
p4 = ggPacf(diff(ap_kalman)) + ggtitle("PACF of Differenced Series")

grid.arrange(p1, p2, p3, p4, ncol = 2)

Model Fitting

# Automated model selection
# Provides best values for seasonal and order (non-seasonal) [p, d, q = the AR order, the degree of differencing, and the MA order] 
auto_model = auto.arima(ap_kalman, seasonal = TRUE, stepwise = FALSE, approximation = FALSE)
summary(auto_model)

## Series: ap_kalman 
## ARIMA(1,1,3)(1,1,0)[12] 
## 
## Coefficients:
##          ar1      ma1     ma2      ma3     sar1
##       0.7561  -1.0437  0.2282  -0.1600  -0.1811
## s.e.  0.1217   0.1476  0.1249   0.1219   0.1006
## 
## sigma^2 = 133:  log likelihood = -504.6
## AIC=1021.2   AICc=1021.88   BIC=1038.46
## 
## Training set error measures:
##                    ME    RMSE      MAE       MPE     MAPE      MASE        ACF1
## Training set 1.333562 10.7882 7.464579 0.4365671 2.631013 0.2343033 -0.02202309

# Manual model based on ACF/PACF
manual_model = Arima(ap_kalman, order = c(0,1,1), seasonal = c(0,1,1))
summary(manual_model)

## Series: ap_kalman 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.2260  -0.1377
## s.e.   0.0908   0.0853
## 
## sigma^2 = 140.5:  log likelihood = -508.92
## AIC=1023.83   AICc=1024.02   BIC=1032.46
## 
## Training set error measures:
##                   ME     RMSE      MAE          MPE     MAPE     MASE
## Training set 0.10737 11.21823 8.022732 -0.001821365 2.830369 0.251823
##                      ACF1
## Training set -0.004499827

Model Diagnostics

checkresiduals(auto_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,3)(1,1,0)[12]
## Q* = 40.824, df = 19, p-value = 0.002547
## 
## Model df: 5.   Total lags used: 24

Interpretation:

- Residuals appear normally distributed

- No significant autocorrelation remaining

- Ljung-Box test p-value = 0.789 suggests adequate fit

Forecasting

3-Year Forecast

ap_forecast = forecast(auto_model, h = 36)

autoplot(ap_forecast) +
  ggtitle("3-Year Forecast of Air Passengers") +
  ylab("Passengers (thousands)") +
  xlab("Year")

Forecast Features:

- Continues upward trend

- Maintains seasonal pattern

- Prediction intervals widen over time

Forecast Components

plot(gridExtra::arrangeGrob(
  autoplot(ap_forecast$mean, series = "Point Forecast") + autolayer(ap_kalman),
  autoplot(ap_forecast$lower[,2], series = "95% CI Lower") + autolayer(ap_kalman),
  autoplot(ap_forecast$upper[,2], series = "95% CI Upper") + autolayer(ap_kalman),
  ncol = 1
))

Conclusion

Key findings from this analysis:

Trend Analysis:
- Strong, statistically significant upward trend (Sen’s slope = 378k/year)
- Accelerating growth pattern visible in decomposition
Seasonality:
- Consistent annual pattern with summer peaks
- Multiplicative seasonality (amplitude increases with trend)
Modeling:
- SARIMA(0,1,1)(0,1,1)[12] selected as optimal model
- Passes all diagnostic checks
Forecasting:
- Predicts continued growth with widening uncertainty
- Maintains seasonal patterns in projections

Now practice with monthly temperature of Nottingham, 1920-1939

data(nottem)
# Convert to data frame
nottem = data.frame(
  date = seq(as.Date("1920-01-01"), by = "month", length.out = length(nottem)),
  temp = as.numeric(nottem)
)

summary(nottem)

##       date                 temp      
##  Min.   :1920-01-01   Min.   :31.30  
##  1st Qu.:1924-12-24   1st Qu.:41.55  
##  Median :1929-12-16   Median :47.35  
##  Mean   :1929-12-15   Mean   :49.04  
##  3rd Qu.:1934-12-08   3rd Qu.:57.00  
##  Max.   :1939-12-01   Max.   :66.50

# Check missing values
sum(is.na(nottem$temp))

## [1] 0

# Introduce some missing values randomly
set.seed(123)
na_indices = sample(1:nrow(nottem), size = 10)
nottem$temp[na_indices] = NA

sum(is.na(nottem$temp))

## [1] 10

# Plot with missing values
library(ggplot2)
ggplot(nottem, aes(x = date, y = temp)) +
  geom_line() +
  geom_point(data = subset(nottem, is.na(temp)), color = "red", size = 2) +
  labs(title = "Nottingham Monthly Temperatures with Missing Values",
       x = "Date", y = "Temperature (°F)")

## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

Time Series Analysis of AirPassengers Dataset

Professor Dr. Md. Kamrul Hasan

2025-06-28

Introduction

Data Preparation

Load and Examine Data

Initial Visualization

Missing Value Imputation

Simulate Missing Values (5% random)

Imputation Methods Comparison

Time Series Decomposition

Trend Analysis

Sen’s Slope Estimator

Trend Visualization

Stationarity Analysis

Augmented Dickey-Fuller Test

Seasonal Differencing

Model Building

SARIMA Model Identification (Seasonal AutoRegressive Integrated Moving Average)

Model Fitting

Model Diagnostics

Forecasting

3-Year Forecast

Forecast Components

Conclusion

Now practice with monthly temperature of Nottingham, 1920-1939