Time Series Regression

Import library

library(tidyverse)
library(gridExtra)
library(readxl)
library(lubridate)
library(hms)
library(forecast)
library(fastDummies)
library(imputeTS)
library(scales)
library(kableExtra)
library(modeltime)
library(Metrics)

Load dataset

df_raw <- read_excel("Data Awal dan Hasil Imputasi rev.xlsx", sheet = 2) %>% 
          mutate(Jam=as_hms(Jam), week = rep(1:51,each=336)[1:16992],
                 time=seq(ymd_hm("2018-1-1 0:00"), ymd_hm("2018-12-20 23:30"), by = 30*60),
                 month=month(time), wday=wday(time),day=day(time), 
                 hour=as_hms(time), yday=yday(time),
                 wtime=rep(seq(ymd_hm("2018-1-1 0:00"), 
                               ymd_hm("2018-1-7 23:30"), 
                               by = 30*60),51)[1:16992]) 
head(df_raw)

var_pollutant <- colnames(df_raw)[6:15]
var_pollutant

 [1] "SUF1_CO"   "SUF7_CO"   "SUF1_NO2"  "SUF7_NO2"  "SUF1_O3"   "SUF7_O3"   "SUF1_PM10" "SUF7_PM10" "SUF1_SO2"  "SUF7_SO2" 

Imputation

df <- list()
df[["na"]] <- data.frame(is.na(df_raw))
df[["impute_1"]] <- df_raw %>% 
                    group_by(Hari, J) %>% 
                    mutate_if(is.double,~replace_na(.,median(., na.rm = TRUE))) %>% 
                    ungroup()

`mutate_if()` ignored the following grouping variables:
• Columns `Hari`, `J`

df[["impute_2"]] <- df_raw %>% 
                    group_by(Hari, J) %>% 
                    mutate_if(is.double,~replace_na(.,mean(., na.rm = TRUE))) %>% 
                    ungroup()

`mutate_if()` ignored the following grouping variables:
• Columns `Hari`, `J`

df[["impute_3"]] <- df_raw %>% 
                    group_by(Hari, J) %>% 
                    mutate_if(is.double,~na_interpolation(., option="linear")) %>% 
                    ungroup()

`mutate_if()` ignored the following grouping variables:
• Columns `Hari`, `J`

Modelling

Dummy variables

df[["dummy"]] <- dummy_cols(df_raw, select_columns = c("J","Hari")) %>% 
                 mutate(t=seq_along(time)) %>% 
                 select(starts_with("J_"), starts_with("Hari_"),t)

for (k in var_pollutant){
  name <- paste0(k,"_imp_median")
  df[[name]] <- bind_cols(df[["impute_1"]], df[["dummy"]]) %>%
                select(time,starts_with(k), starts_with("J_"), starts_with("Hari_"),t)
  name <- paste0(k,"_imp_mean")
  df[[name]] <- bind_cols(df[["impute_2"]], df[["dummy"]]) %>% 
                select(time,starts_with(k), starts_with("J_"), starts_with("Hari_"),t)
  name <- paste0(k,"_imp_interpolation")
  df[[name]] <- bind_cols(df[["impute_3"]], df[["dummy"]]) %>% 
                select(time,starts_with(k), starts_with("J_"), starts_with("Hari_"),t)
}

ind <- df_raw %>% make_ts_splits(.length_test = 48)
ind

split [16944|48|16992]

Model Fitting and Forecasting

applyForecast <- function(data){
  ind <- data %>% make_ts_splits(.length_test = 48)
  train <- data[ind$idx_train,]
  test  <- data[ind$idx_test,]
  yy <- colnames(data)[2]
  f <- paste0(yy,"~.-1-time-Hari_1")
  m <- lm(f, data = train)
  r.squared <- summary(m)$r.squared
  pred <- m %>% predict(test)
  RMSE <- rmse(test[yy] %>% pull(),pred)
  MAPE <- mape(test[yy] %>% pull(),pred)
  sMAPE <- smape(test[yy] %>% pull(),pred)
  p_train <- ggplot(train %>% mutate(predicted=m$fitted.values)) + 
             geom_line(aes_(x=~time,y=as.name(yy))) +  
             geom_line(aes(x=time,y=predicted), color="red") +
             xlab("") +
             scale_x_datetime(date_breaks = "4 hour", date_labels = "%b %e %H:%M")
  p_test  <- ggplot(test %>% mutate(predicted=pred)) + 
             geom_line(aes_(x=~time,y=as.name(yy))) +  
             geom_line(aes(x=time,y=predicted), color="red") +
             xlab("") +
             scale_x_datetime(date_breaks = "4 hour", date_labels = "%b %e %H:%M")
  results <- list(prediction=pred, 
                  r.squared=r.squared,
                  RMSE=RMSE,
                  MAPE=MAPE,
                  sMAPE=sMAPE,
                  plot_train=p_train,
                  plot_test=p_test)
}

df_pollutant <- names(df)[grepl("imp_",names(df))]
res <- list()
for (name in df_pollutant){res[[name]] <- applyForecast(df[[name]])}

Forecast Comparison

dd <- as_tibble(res[[1]][2:5])
for (name in df_pollutant[-1]){
  dd <- dd %>% add_row(as_tibble(res[[name]][2:5]))
}
results <- dd %>% mutate(model=df_pollutant)

CO

results[1:6,]

grid.arrange(arrangeGrob(res[[2]]$plot_test, top = "TSR with Mean Imputation"),
             arrangeGrob(res[[4]]$plot_test, top = "TSR with Median Imputation"),
             ncol=1)

NO2

results[7:12,]

grid.arrange(arrangeGrob(res[[8]]$plot_test, top = "TSR with Mean Imputation"),
             arrangeGrob(res[[10]]$plot_test, top = "TSR with Median Imputation"),
             ncol=1)

O3

results[13:18,]

grid.arrange(arrangeGrob(res[[14]]$plot_test, top = "TSR with Mean Imputation"),
             arrangeGrob(res[[18]]$plot_test, top = "TSR with Linear Interpolation Imputation"),
             ncol=1)

PM10

results[19:24,]

grid.arrange(arrangeGrob(res[[20]]$plot_test, top = "TSR with Mean Imputation"),
             arrangeGrob(res[[22]]$plot_test, top = "TSR with Median Imputation"),
             ncol=1)

SO2

results[25:30,]

grid.arrange(arrangeGrob(res[[26]]$plot_test, top = "TSR with Mean Imputation"),
             arrangeGrob(res[[28]]$plot_test, top = "TSR with Median Imputation"),
             ncol=1)