Short-term probabilistic forecasting of GB national net demand

for (name in names(list_of_forms)){
  if(file.exists(
    paste0("~/RProjects/reactforecasting/data/output/R_objects/",
           "ShortTerm_NetDemand/qgam_update-",name,".rds"))){
    print(paste0("Model '", name,"' already fitted."))
    
  }else{
    start_time <- Sys.time()
    print(paste0("Model: ", name))
    selected_list_of_forms <- list_of_forms[[name]]
    selected_res_sq_form <- list_of_res_sq_forms[[name]]
    selected_qr_form <- list_of_qr_forms[[name]]
    qgam <- qreg_gam_ma(data = fit_data,
                        formula_list = selected_list_of_forms,
                        response = "net_demand",
                        model_r_sq = T,
                        formula_r_sq = selected_res_sq_form,
                        kfold = "kfold",
                        exclude_train = 'missing_data',
                        formula_qr = selected_qr_form,
                        weights = weights,
                        quantiles = c(0.025, seq(0.05, 0.95, 0.05),
                                      0.975)
                        )
  
    qgam_update <- copy(qgam)
  
    saveRDS(qgam, 
            paste0("~/RProjects/reactforecasting/data/output/R_objects/",
            "ShortTerm_NetDemand/qgam-",name,".rds"))
    
    rm(qgam)
    test_issues <- fit_data[kfold=="Test",(unique(issueTime))]
    weekly_test_issues <- (split(split(test_issues, test_issues),
                                 ceiling(seq_along(test_issues)/(14*48))))
    for(week in 1:(length(weekly_test_issues)-1)){
      week_index = which(fit_data$issueTime%in%weekly_test_issues[[week]]
                         &fit_data$kfold=="Test" )
      next_week_index = which(fit_data$issueTime%in%weekly_test_issues[[week+1]]
                              &fit_data$kfold=="Test" )
      
      previous_day <- fit_data[week_index,min(issueTime)] - days(1)
      
      last_day <- fit_data[next_week_index,max(issueTime)]
      
      fit_data_week <- fit_data[issueTime>=previous_day &
                                  issueTime <= last_day]
      
      week_index_sub = which(fit_data_week$issueTime%in%weekly_test_issues[[week]]
                         &fit_data_week$kfold=="Test" )
      next_week_index_sub = which(fit_data_week$issueTime%in%weekly_test_issues[[week+1]]
                              &fit_data_week$kfold=="Test" )
      
      print(paste0("Updating forecast for 4-week block ",
                   week," of ",length(weekly_test_issues)-1))
    
      ## Update "Test" model with previous block's data
      qgam_update <- qreg_gam_ma.update(qgam_update, data=fit_data_week, 
                                        update_index = week_index_sub)
    
      ## Make new predictions for current block
      new_preds <- predict(qgam_update,data = fit_data_week, 
                           predict_index = next_week_index_sub, predict_quantiles = T)
    
      ## Update test predictions for current block
      qgam_update$mqr_pred[next_week_index,] <- new_preds$mqr_pred
      qgam_update$models$gam_pred$gam_pred[next_week_index] <- new_preds$gam_pred
      qgam_update$models$gam_pred$gam_sd[next_week_index] <- new_preds$gam_sd
    
      rm(new_preds)
      if(name=="model"){
        print("All lead times.")
        fit_data_all_leads <- readRDS("~/RProjects/reactforecasting/data/output/R_objects/ShortTerm_NetDemand/fit_data_all_leads.rds")
        fit_data_all_week <- fit_data_all_leads[issueTime>=previous_day &
                                                  issueTime <= last_day]
        next_week_index_all_sub = which(fit_data_all_week$issueTime 
                                        %in% weekly_test_issues[[week+1]]
                                        & fit_data_all_week$kfold=="Test" )
        next_week_index_all = which(fit_data_all_leads$issueTime 
                                    %in% weekly_test_issues[[week+1]]
                                    & fit_data_all_leads$kfold=="Test" )
        all_new_preds <- predict(qgam_update,
                                 data = fit_data_all_week, 
                                 predict_index = next_week_index_all_sub, 
                                 predict_quantiles = T)
        rm(fit_data_all_leads)
        
        pred_all_leads <- readRDS("~/RProjects/reactforecasting/data/output/R_objects/ShortTerm_NetDemand/pred_all_leads.rds")
        pred_all_leads$mqr_pred[next_week_index_all] <- all_new_preds$mqr_pred
        pred_all_leads$gam_pred[next_week_index_all] <- all_new_preds$gam_pred
        pred_all_leads$gam_sd[next_week_index_all] <- all_new_preds$gam_sd
        rm(all_new_preds)
        saveRDS(pred_all_leads, 
                "~/RProjects/reactforecasting/data/output/R_objects/ShortTerm_NetDemand/pred_all_leads.rds")
        rm(pred_all_leads)
      }
    }
    
    saveRDS(qgam_update, 
            paste0("~/RProjects/reactforecasting/data/output/R_objects/",
            "ShortTerm_NetDemand/qgam_update-",name,".rds"))
    
    rm(qgam_update)
    rm(test_issues)
    rm(weekly_test_issues)
    rm(week_data)
    rm(next_week_data)
    rm(selected_list_of_forms)
    rm(selected_res_sq_form)
    print(paste0("Model fitting time for ", name,": ", format(Sys.time() - start_time)))
  }
}