Forecasting_Final
Yesh Onipede
2024-04-24
Downloading the data
rm(list = ls())
setwd('/Users/yeshimonipede/Desktop/BC_Spring2024')
train_data<-read.csv("dengue_labels_train.csv")
test_data<-read.csv("dengue_features_test.csv")
ex_regressors<-read.csv("dengue_features_train.csv")
submission_file<-read.csv("submission_format.csv")
#Different train and test sets for each location
train_sj <- train_data[train_data$city == "sj", ]
train_iq<-train_data[train_data$city == "iq", ]
test_sj<-test_data[test_data$city=="sj",]
test_iq<-test_data[test_data$city == "iq", ]
ex_regressors_sj<-ex_regressors[ex_regressors$city=="sj",]
ex_regressors_iq<-ex_regressors[ex_regressors$city=="iq",]
library(ggplot2)
library(tidyr)
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericDealing with NAs - fills the missing values in each train and test set with the most recent values along each column
train_sj_filled <- na.locf(train_sj)
train_iq_filled <- na.locf(train_iq)
test_sj_filled <- na.locf(test_sj)
test_iq_filled <- na.locf(test_iq)
ex_regressors_sj_filled <- na.locf(ex_regressors_sj)
ex_regressors_iq_filled <- na.locf(ex_regressors_iq)Turning sj datasets in ts
#This is the time series data for the entire sj train dataset provided
ts_train_sj_filled <- ts(train_sj_filled[, -c(1:2)],
start = c(1990, 18),
end = c(2010, 17),
frequency = 53)
sj_train_univariate <- ts(train_sj_filled[, 4],
start = c(1990, 18),
end = c(2008, 17),
frequency = 53)
sj_test <- window(ts_train_sj_filled,
start=c(2008,18),
frequency=53)
sj_test_univariate <- ts(sj_test[,2],
start=c(2008,18),
frequency=53)Plotting the data
# Plot the sj_train_univariate time series
plot(sj_train_univariate, main = "Train Data - San Juan (Univariate)", xlab = "Time", ylab = "Values")
STL the training data
# Decompose the time series using STL
train_sj_stl <- stl(sj_train_univariate, s.window = "periodic")
# Plot the decomposition
plot(train_sj_stl)
Arima model - all train data
# Fit ARIMA model to the univariate time series
sj_arima_model <- forecast::auto.arima(sj_train_univariate)
# Make forecasts for the test data using the fitted ARIMA model
sj_arima_forecast <- forecast::forecast(sj_arima_model, h = 106)
sj_arima_forecast_submit <- forecast::forecast(sj_arima_model, h = 275)
# Print the ARIMA model summary
print(summary(sj_arima_model))## Series: sj_train_univariate
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.7127 -0.5947
## s.e. 0.0936 0.1064
##
## sigma^2 = 178.1: log likelihood = -3820.56
## AIC=7647.12 AICc=7647.15 BIC=7661.7
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.0211412 13.32308 7.977775 NaN Inf 0.2190721 0.0011972accuracy(sj_arima_forecast,sj_test_univariate)## ME RMSE MAE MPE MAPE MASE
## Training set 0.0211412 13.32308 7.977775 NaN Inf 0.2190721
## Test set 15.2928976 39.56351 25.940882 -7.338393 53.04048 0.7123445
## ACF1 Theil's U
## Training set 0.0011972 NA
## Test set 0.9290522 2.141422Plotting the residuals: The residuals do not necessarily exhibit white noise well which leads me to believe that an ARIMA model is not the best model for this dataset
# Obtain the residuals from the fitted ARIMA model
sj_arima_residuals <- residuals(sj_arima_model)
# Create a time series object for the residuals
sj_arima_residuals_ts <- ts(sj_arima_residuals, start = start(sj_train_univariate))
# Plot the residuals over time
plot(sj_arima_residuals_ts, main = "Residuals of ARIMA Model", xlab = "Time", ylab = "Residuals")
Neural Net Model - all train data
# Fit neural network model using nnetar()
sj_nnet_model <- forecast::nnetar(sj_train_univariate)
# Make forecasts for the test data using the fitted neural network model
sj_nnet_forecast_submit <- forecast::forecast(sj_nnet_model, h = 275)
sj_nnet_forecast <- forecast::forecast(sj_nnet_model, h = 106)
# Print the neural network model summaryArima with Exogenous Regressors
# Fit ARIMA model to the univariate time series
ex_regs <- ex_regressors_sj_filled[, c("station_avg_temp_c","reanalysis_specific_humidity_g_per_kg","reanalysis_dew_point_temp_k","station_min_temp_c")]
names(ex_regs) <- c("reg1", "reg2", "reg3", "reg4")
# Create separate vectors for each column
reg1 <- ex_regs[, "reg1"]
reg2 <- ex_regs[, "reg2"]
reg3 <- ex_regs[, "reg3"]
reg4 <- ex_regs[, "reg4"]
# Convert each column to a time series
reg1_ts <- ts(reg1, start = c(1990, 18), end = c(2008, 17), frequency = 53)
reg2_ts <- ts(reg2, start = c(1990, 18), end = c(2008, 17), frequency = 53)
reg3_ts <- ts(reg3, start = c(1990, 18), end = c(2008, 17), frequency = 53)
reg4_ts <- ts(reg4, start = c(1990, 18), end = c(2008, 17), frequency = 53)
# Fit ARIMA model to the univariate time series with exogenous regressors
sj_arima_model_ex <- forecast::auto.arima(sj_train_univariate, xreg = cbind(reg1_ts, reg2_ts, reg3_ts, reg4_ts))
# Make forecasts for the test data using the fitted ARIMA model with exogenous regressors
sj_ex_submit_forecast <- forecast::forecast(sj_arima_model_ex, xreg = cbind(reg1_ts, reg2_ts, reg3_ts, reg4_ts), h = 275)
sj_ex_forecast <- forecast::forecast(sj_arima_model_ex, xreg = cbind(reg1_ts, reg2_ts, reg3_ts, reg4_ts), h = 106)Plotting forecasts - sj
autoplot(sj_train_univariate,ylab = "Total Cases")+
autolayer(sj_nnet_forecast,series="Neural Net")+
autolayer(sj_arima_forecast,series="ARIMA model",PI=FALSE)+
autolayer(sj_ex_forecast,series="Regression with Exogenous Regressor",PI=FALSE)
# Select data from the year 2000 onwards
sj_train_univariate_zoomed <- window(sj_train_univariate, start = c(2005, 1))
# Create autoplot for the zoomed-in data with y-axis labeled as "Total Cases"
autoplot(sj_train_univariate_zoomed, ylab = "Total Cases") +
autolayer(sj_nnet_forecast, series = "Neural Net") +
autolayer(sj_arima_forecast, series = "ARIMA model", PI = FALSE) +
autolayer(sj_ex_forecast, series = "Regression with Exogenous Regressor", PI = FALSE)
Accuracy against test set
# ARIMA model accuracy
arima_accuracy <- accuracy(sj_arima_forecast, sj_test_univariate)
print("ARIMA Model Accuracy:")## [1] "ARIMA Model Accuracy:"print(arima_accuracy)## ME RMSE MAE MPE MAPE MASE
## Training set 0.0211412 13.32308 7.977775 NaN Inf 0.2190721
## Test set 15.2928976 39.56351 25.940882 -7.338393 53.04048 0.7123445
## ACF1 Theil's U
## Training set 0.0011972 NA
## Test set 0.9290522 2.141422# Neural Network model accuracy
nnet_accuracy <- accuracy(sj_nnet_forecast, sj_test_univariate)
print("Neural Network Model Accuracy:")## [1] "Neural Network Model Accuracy:"print(nnet_accuracy)## ME RMSE MAE MPE MAPE MASE
## Training set 0.0287018 7.120509 5.150439 -Inf Inf 0.1414326
## Test set 6.3263355 56.722356 41.916582 -77.78106 135.2131 1.1510420
## ACF1 Theil's U
## Training set -0.007281439 NA
## Test set 0.967170502 7.517314exreg_accuracy <- accuracy(sj_ex_submit_forecast , sj_test_univariate)
print("Exogenous Regressor Model Accuracy:")## [1] "Exogenous Regressor Model Accuracy:"print(exreg_accuracy)## ME RMSE MAE MPE MAPE MASE
## Training set 0.01967364 13.23152 7.990913 -Inf Inf 0.2194329
## Test set 19.80469625 40.76393 25.679137 9.718971 44.75959 0.7051569
## ACF1 Theil's U
## Training set -0.000640189 NA
## Test set 0.924144874 1.834487Based on the test set metrics, the Exogenous Regressor Model appears to have performed the best on the test set, as it has the lowest RMSE, MAE, and MAPE among the three models.
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
IQ DATA
ts_train_iq_filled <- ts(train_iq_filled[, -c(1:2)],
start = c(1990, 18),
end = c(2010, 17),
frequency = 53)
iq_train_univariate <- ts(train_iq_filled[, 4],
start = c(1990, 18),
end = c(2008, 17),
frequency = 53)
iq_test <- window(ts_train_iq_filled,
start=c(2008,18),
frequency=53)
iq_test_univariate <- ts(iq_test[,2],
start=c(2008,18),
frequency=53)# Plot the iq_train_univariate time series
plot(iq_train_univariate, main = "Train Data - Iquitos (Univariate)", xlab = "Time", ylab = "Values")
# Decompose the time series using STL
train_sj_stl <- stl(sj_train_univariate, s.window = "periodic")
# Plot the decomposition
plot(train_sj_stl)
ARIMA Model
# Fit ARIMA model to the univariate time series
iq_arima_model <- forecast::auto.arima(iq_train_univariate)
# Make forecasts for the test data using the fitted ARIMA model
iq_arima_forecast <- forecast::forecast(iq_arima_model, h = 106)
iq_arima_forecast_submit <- forecast::forecast(iq_arima_model, h = 275)
# Print the ARIMA model summary
print(summary(iq_arima_model))## Series: iq_train_univariate
## ARIMA(3,0,2) with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 mean
## 0.8675 -0.6197 0.5204 -0.1849 0.4289 7.8004
## s.e. 0.1003 0.1258 0.0645 0.1107 0.1002 1.2257
##
## sigma^2 = 50.62: log likelihood = -3223.16
## AIC=6460.32 AICc=6460.44 BIC=6494.34
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.008496431 7.092705 3.917964 -Inf Inf 0.4237289 -0.005019387accuracy(iq_arima_forecast,iq_test_univariate)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.008496431 7.092705 3.917964 -Inf Inf 0.4237289 -0.005019387
## Test set -3.944447005 8.090438 6.647624 -Inf Inf 0.7189424 0.819055650
## Theil's U
## Training set NA
## Test set 0Neural Net Model
# Fit neural network model using nnetar()
iq_nnet_model <- forecast::nnetar(iq_train_univariate)
# Make forecasts for the test data using the fitted neural network model
iq_nnet_forecast_submit <- forecast::forecast(iq_nnet_model, h = 275)
iq_nnet_forecast <- forecast::forecast(iq_nnet_model, h = 106)Exogenous Regressor Model
# Fit ARIMA model to the univariate time series
ex_regs <- ex_regressors_iq_filled[, c("station_avg_temp_c","reanalysis_specific_humidity_g_per_kg","reanalysis_dew_point_temp_k","station_min_temp_c")]
names(ex_regs) <- c("reg1", "reg2", "reg3", "reg4")
# Create separate vectors for each column
reg1 <- ex_regs[, "reg1"]
reg2 <- ex_regs[, "reg2"]
reg3 <- ex_regs[, "reg3"]
reg4 <- ex_regs[, "reg4"]
# Convert each column to a time series
reg1_ts <- ts(reg1, start = c(1990, 18), end = c(2008, 17), frequency = 53)
reg2_ts <- ts(reg2, start = c(1990, 18), end = c(2008, 17), frequency = 53)
reg3_ts <- ts(reg3, start = c(1990, 18), end = c(2008, 17), frequency = 53)
reg4_ts <- ts(reg4, start = c(1990, 18), end = c(2008, 17), frequency = 53)
# Fit ARIMA model to the univariate time series with exogenous regressors
iq_arima_model_ex <- forecast::auto.arima(iq_train_univariate, xreg = cbind(reg1_ts, reg2_ts, reg3_ts, reg4_ts))
# Make forecasts for the test data using the fitted ARIMA model with exogenous regressors
iq_ex_submit_forecast <- forecast::forecast(iq_arima_model_ex, xreg = cbind(reg1_ts, reg2_ts, reg3_ts, reg4_ts), h = 275)
iq_ex_forecast <- forecast::forecast(iq_arima_model_ex, xreg = cbind(reg1_ts, reg2_ts, reg3_ts, reg4_ts), h = 106)# Plotting forecasts for the IQ data
autoplot(iq_train_univariate, ylab = "Total Cases") +
autolayer(iq_nnet_forecast, series = "Neural Net") +
autolayer(iq_arima_forecast, series = "ARIMA model", PI = FALSE) +
autolayer(iq_ex_forecast, series = "Regression with Exogenous Regressor", PI = FALSE)
# Select data from the year 2000 onwards
iq_train_univariate_zoomed <- window(iq_train_univariate, start = c(2005, 1))
# Create autoplot for the zoomed-in data with y-axis labeled as "Total Cases"
autoplot(iq_train_univariate_zoomed, ylab = "Total Cases") +
autolayer(iq_nnet_forecast, series = "Neural Net") +
autolayer(iq_arima_forecast, series = "ARIMA model", PI = FALSE, size=1.2) +
autolayer(iq_ex_forecast, series = "Regression with Exogenous Regressor", PI = FALSE)
# ARIMA model accuracy for IQ data
iq_arima_accuracy <- accuracy(iq_arima_forecast, iq_test_univariate)
print("ARIMA Model Accuracy:")## [1] "ARIMA Model Accuracy:"print(iq_arima_accuracy)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.008496431 7.092705 3.917964 -Inf Inf 0.4237289 -0.005019387
## Test set -3.944447005 8.090438 6.647624 -Inf Inf 0.7189424 0.819055650
## Theil's U
## Training set NA
## Test set 0# Neural Network model accuracy for IQ data
iq_nnet_accuracy <- accuracy(iq_nnet_forecast, iq_test_univariate)
print("Neural Network Model Accuracy:")## [1] "Neural Network Model Accuracy:"print(iq_nnet_accuracy)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.04260568 3.665107 2.584685 NaN Inf 0.2795345 -0.01587416
## Test set -7.81981788 12.495573 10.408844 -Inf Inf 1.1257194 0.88254338
## Theil's U
## Training set NA
## Test set 0# Exogenous Regressor model accuracy for IQ data
iq_exreg_accuracy <- accuracy(iq_ex_forecast, iq_test_univariate)
print("Exogenous Regressor Model Accuracy:")## [1] "Exogenous Regressor Model Accuracy:"print(iq_exreg_accuracy)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.01269759 7.084116 3.931561 NaN Inf 0.4251994 -0.004751527
## Test set -4.18112717 8.256624 6.777280 -Inf Inf 0.7329648 0.817036683
## Theil's U
## Training set NA
## Test set 0Based on the accuracy metrics for the test set, the ARIMA Model appears to have performed the best among the three models for the IQ dataset. It has the lowest RMSE, MAE, and MAPE values compared to the other models.
Plotting residuals of ARIMA
# Obtain the residuals from the fitted ARIMA model
iq_arima_residuals <- residuals(iq_arima_model)
# Create a time series object for the residuals
iq_arima_residuals_ts <- ts(iq_arima_residuals, start = start(iq_train_univariate))
# Plot the residuals over time
plot(iq_arima_residuals_ts, main = "Residuals of ARIMA Model", xlab = "Time", ylab = "Residuals")