Time Series Imputation

Libraries

library(tidyverse)
library(imputeTS)
library(patchwork)
library(seastests)
library(trend)
library(forecast)

Dataset

link <- "https://gist.githubusercontent.com/novrisuhermi/cc55e030a8d73a6b286976bd8d48dec2/raw/e610da2a0043e260d78692e2f82e1c44ea936ee4/ElectricProduction.csv"
df <- read.csv(link)
df$DATE <- as.Date(df$DATE, format = "%m/%d/%Y")

p <- ggplot(df, aes(x = DATE, y = IPG2211A2N)) +
  geom_line(colour = "steelblue4", linewidth = 0.5) +
  scale_x_date(
    breaks = seq(as.Date("1985-01-01"), as.Date("2015-01-01"), by = "5 years"),
    date_labels = "%Y"
  ) +
  labs(
    title    = "Monthly Electric Production",
    subtitle = "January 1985 - January 2018",
    x        = "Year",
    y        = "Production index"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title       = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

p

Seasonality and Trend

y <- ts(df$IPG2211A2N, start = c(1985, 1), frequency = 12)

Seasonal Plot

p_season <- ggseasonplot(y, continuous = TRUE) +
  labs(title = "Seasonal Plot", x = "Month", y = "Production index", colour = "Year") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

p_subseries <- ggsubseriesplot(y) +
  labs(title = "Seasonal Subseries Plot", subtitle = "Blue lines = monthly means",
       x = "Month", y = "Production index") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

p_season | p_subseries

Seasonality Tests

# H0: no seasonality. By default the tests use the first-differenced series,
# which removes the trend.
kw(y)                              # Kruskal-Wallis test

Test used:  Kruskall Wallis 
 
Test statistic:  345.35 
P-value:  0

fried(y)                           # Friedman test

Test used:  Friedman rank 
 
Test statistic:  323.19 
P-value:  0

qs(y)                              # QS test (autocorrelation at seasonal lags)

Test used:  QS 
 
Test statistic:  610.31 
P-value:  0

isSeasonal(y, test = "combined")   # Webel-Ollech combined test: TRUE = seasonal

[1] TRUE

Trend Tests

# H0: no monotonic trend
smk.test(y)   # Seasonal Mann-Kendall (accounts for seasonality)

    Seasonal Mann-Kendall trend test (Hirsch-Slack test)

data:  y
z = 24.471, p-value < 2.2e-16
alternative hypothesis: true S is not equal to 0
sample estimates:
    S  varS 
 5493 50369 

mk.test(y)    # Ordinary Mann-Kendall, for comparison

    Mann-Kendall trend test

data:  y
z = 19.125, n = 397, p-value < 2.2e-16
alternative hypothesis: true S is not equal to 0
sample estimates:
           S         varS          tau 
5.052400e+04 6.978466e+06 6.427499e-01 

# Seasonal Sen's slope = change per year; needs complete years, so Jan 2018 is dropped
sea.sens.slope(window(y, end = c(2017, 12)))

[1] 1.392249

Simulating Missing Values

inject_gaps <- function(x, prop = 0.15, gap_len = 1:40, gap_prob = 1 / gap_len,
                        cover = list(single = 1, short = 2:5, long = 6:12),
                        min_sep = 1, buffer = 1) {
  n <- length(x)
  m <- round(prop * n)

  # 1. Gap lengths: one per class in `cover`, the rest random, summing to m
  L <- vapply(cover, function(r) r[sample.int(length(r), 1)], numeric(1))
  if (sum(L) > m) stop("prop is too small for the gap classes in `cover`.")
  L <- c(L, gap_len[sample.int(length(gap_len), m, replace = TRUE, prob = gap_prob)])
  k <- which(cumsum(L) >= m)[1]
  L <- L[seq_len(k)]
  L[k] <- m - sum(L[-k])
  L <- unname(L[sample.int(k)])

  # 2. Observed points that are free to be placed anywhere
  free <- n - m - 2 * buffer - (k - 1) * min_sep
  if (free < 0) stop("Gaps do not fit: reduce prop, min_sep or buffer.")

  # 3. Spread the free points randomly over the k + 1 spaces (stars and bars)
  bars  <- sort(sample.int(free + k, k))
  extra <- diff(c(0, bars, free + k + 1)) - 1

  # 4. Gap start positions, then set the gaps to NA
  space <- c(buffer, rep(min_sep, k - 1)) + extra[seq_len(k)]
  start <- cumsum(space + c(0, L[-k])) + 1
  x[sequence(L, from = start)] <- NA
  x
}

set.seed(2026)                                   
df_missing <- df %>%
  mutate(IPG2211A2N = inject_gaps(IPG2211A2N, prop = 0.20))

df_missing

# Left: series with missing regions (x = date)
p_series <- ggplot_na_distribution(df_missing$IPG2211A2N,
                                   x_axis_labels = df_missing$DATE,
                                   xlab = "Year",
                                   ylab = "Production index") +
  scale_x_date(
    breaks = seq(as.Date("1985-01-01"), as.Date("2015-01-01"), by = "5 years"),
    date_labels = "%Y"
  )

# Right: share of NA per year (x = observation index, relabelled as years)
p_interval <- suppressMessages(
  ggplot_na_distribution2(df_missing$IPG2211A2N, interval_size = 12, xlab = "Year") +
    scale_x_continuous(breaks = seq(0, 360, by = 60),     # index 0, 60, ... = Jan 1985, Jan 1990, ...
                       labels = seq(1985, 2015, by = 5),
                       expand = c(0, 0))
)

p_series | p_interval

Imputation

y_na <- ts(df_missing$IPG2211A2N, start = c(1985, 1), frequency = 12)

quiet_structts <- function(expr) {
  withCallingHandlers(expr, warning = function(w) {
    if (grepl("possible convergence problem", conditionMessage(w))) invokeRestart("muffleWarning")
  })
}

set.seed(2026)   # na_random() and the "random" algorithms draw random values
imputations <- list(
  # Mean-based
  mean      = na_mean(y_na, option = "mean"),
  median    = na_mean(y_na, option = "median"),
  mode      = na_mean(y_na, option = "mode"),
  harmonic  = na_mean(y_na, option = "harmonic"),
  geometric = na_mean(y_na, option = "geometric"),

  # Random values between the observed minimum and maximum
  random = na_random(y_na),

  # Last observation carried forward / next observation carried backward
  locf = na_locf(y_na, option = "locf"),
  nocb = na_locf(y_na, option = "nocb"),

  # Interpolation
  linear = na_interpolation(y_na, option = "linear"),
  spline = na_interpolation(y_na, option = "spline"),
  stine  = na_interpolation(y_na, option = "stine"),

  # Weighted moving average (k = 4 observations on each side)
  ma_simple      = na_ma(y_na, k = 4, weighting = "simple"),
  ma_linear      = na_ma(y_na, k = 4, weighting = "linear"),
  ma_exponential = na_ma(y_na, k = 4, weighting = "exponential"),

  # Kalman smoothing
  kalman_structts = quiet_structts(na_kalman(y_na, model = "StructTS")),
  kalman_arima    = na_kalman(y_na, model = "auto.arima"),

  # Seasonal decomposition: impute the deseasonalised series, then add the season back
  seadec_interpolation = na_seadec(y_na, algorithm = "interpolation"),
  seadec_locf          = na_seadec(y_na, algorithm = "locf"),
  seadec_mean          = na_seadec(y_na, algorithm = "mean"),
  seadec_random        = na_seadec(y_na, algorithm = "random"),
  seadec_kalman        = quiet_structts(na_seadec(y_na, algorithm = "kalman")),
  seadec_ma            = na_seadec(y_na, algorithm = "ma"),

  # Seasonal split: impute each calendar month's subseries separately
  seasplit_interpolation = na_seasplit(y_na, algorithm = "interpolation"),
  seasplit_locf          = na_seasplit(y_na, algorithm = "locf"),
  seasplit_mean          = na_seasplit(y_na, algorithm = "mean"),
  seasplit_random        = na_seasplit(y_na, algorithm = "random"),
  seasplit_kalman        = quiet_structts(na_seasplit(y_na, algorithm = "kalman")),
  seasplit_ma            = na_seasplit(y_na, algorithm = "ma")
)

df_imputed <- df_missing %>%
  select(DATE, observed = IPG2211A2N) %>%
  mutate(actual = df$IPG2211A2N, .after = DATE) %>%
  bind_cols(as_tibble(map(imputations, as.numeric)))

df_imputed

Evaluation

RMSE and MAPE

Accuracy is computed only on the m imputed (originally missing) points \mathcal{M}:

\text{RMSE} = \sqrt{\frac{1}{m}\sum_{t \in \mathcal{M}} (\hat{y}_t - y_t)^2}, \qquad \text{MAPE} = \frac{100}{m}\sum_{t \in \mathcal{M}} \left\lvert \frac{\hat{y}_t - y_t}{y_t} \right\rvert

method_groups <- list(
  "Mean-based"             = c("mean", "median", "mode", "harmonic", "geometric"),
  "Random"                 = "random",
  "LOCF / NOCB"            = c("locf", "nocb"),
  "Interpolation"          = c("linear", "spline", "stine"),
  "Moving Average"         = c("ma_simple", "ma_linear", "ma_exponential"),
  "Kalman Smoothing"       = c("kalman_structts", "kalman_arima"),
  "Seasonal Decomposition" = paste0("seadec_",   c("interpolation", "locf", "mean", "random", "kalman", "ma")),
  "Seasonal Split"         = paste0("seasplit_", c("interpolation", "locf", "mean", "random", "kalman", "ma"))
)

# RMSE and MAPE on the imputed points only
evaluation <- df_imputed %>%
  filter(is.na(observed)) %>%
  pivot_longer(-c(DATE, actual, observed), names_to = "method", values_to = "imputed") %>%
  group_by(method) %>%
  summarise(RMSE = sqrt(mean((imputed - actual)^2)),
            MAPE = 100 * mean(abs(imputed - actual) / actual)) %>%
  left_join(enframe(method_groups, name = "group", value = "method") %>% unnest(method),
            by = "method") %>%
  relocate(group) %>%
  arrange(RMSE)

# Display: ranked table, rounded, all 28 rows on one page (chunk option rows.print)
evaluation %>%
  mutate(Rank = row_number(), across(c(RMSE, MAPE), ~ round(.x, 3))) %>%
  select(Rank, Family = group, Method = method, RMSE, `MAPE (%)` = MAPE)

Imputed Time Series Visualization

plot_imputations <- function(group, ncol = 2) {
  methods <- method_groups[[group]]
  labels  <- evaluation %>%
    mutate(label = sprintf("%s  (RMSE %.2f, MAPE %.2f%%)", method, RMSE, MAPE))

  df_imputed %>%
    select(DATE, actual, observed, all_of(methods)) %>%
    pivot_longer(all_of(methods), names_to = "method", values_to = "imputed") %>%
    mutate(method  = factor(method, levels = methods,
                            labels = labels$label[match(methods, labels$method)]),
           imputed = if_else(is.na(observed), imputed, NA_real_)) %>%
    ggplot(aes(x = DATE)) +
    geom_line(aes(y = actual), colour = "grey70", linewidth = 0.3) +
    geom_line(aes(y = observed), colour = "steelblue4", linewidth = 0.4, na.rm = TRUE) +
    geom_point(aes(y = imputed), colour = "firebrick", size = 0.8, na.rm = TRUE) +
    facet_wrap(~ method, ncol = ncol) +
    scale_x_date(
      breaks = seq(as.Date("1985-01-01"), as.Date("2015-01-01"), by = "5 years"),
      date_labels = "%Y"
    ) +
    labs(title    = group,
         subtitle = "Grey = actual (hidden) values, blue = observed, red = imputed",
         x = "Year", y = "Production index") +
    theme_minimal(base_size = 11) +
    theme(plot.title = element_text(face = "bold"))
}

Mean-based

plot_imputations("Mean-based")

Random

plot_imputations("Random", ncol = 1)

LOCF / NOCB

plot_imputations("LOCF / NOCB")

Interpolation

plot_imputations("Interpolation")

Moving Average

plot_imputations("Moving Average")

Kalman Smoothing

plot_imputations("Kalman Smoothing")

Seasonal Decomposition

plot_imputations("Seasonal Decomposition")

Seasonal Split

plot_imputations("Seasonal Split")