general_model_selection.Rmd

---
title: "Var_TS"
author: "Jaskirat Kaur"
date: "2023-11-30"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, include=FALSE}
library(tidyverse)
library(forecast)
library(sf)
library(imputeTS)
library(data.table)
library(tseries)
```


### Loading data

```{r}
lakes <- fread('/Users/jaskiratkaur/Documents/TimeSeries/data/aggregated_lakes.csv')
lakes
```

## Time series
```{r}
ts_water <- ts(lakes$water_area, frequency = 12) %>% na_interpolation()
ts_evap <- ts(lakes$evap_volume, frequency = 12) %>% na_interpolation()
```


```{r}
tsdisplay(ts_water)
```

### ARIMA

```{r}
arima_water_area <- auto.arima(ts_water, lambda="auto")
summary(arima_water_area)
```
```{r}
checkresiduals(arima_water_area)
```

```{r}
plot(lakes$water_area, lakes$evap_volume, main = "Water Area vs Evaporation volume", xlab = "Water Area", ylab = "Evaporation Volume")
```

```{r}
correlation <- cor(lakes$water_area, lakes$evap_volume)
cat("Correlation between Water Area and Evaporation volume:", correlation, "\n")
```

### ARIMAX (with evaporation volume)

```{r}
arimax_model <- auto.arima(ts_water, xreg = ts_evap, lambda="auto")
summary(arimax_model)
```
```{r}
checkresiduals(arimax_model)
```


### GARCH

```{r}
garch_model <- garch(ts_water)

```

```{r}
AIC(garch_model)
```

### VAR (with evaporation volume)


```{r}
# Fit VAR(10) model
var_model <- VAR(cbind(lakes$water_area, lakes$evap_volume), p = 12, type = "both", season = 12)

```


```{r}
checkresiduals(var_model[,1])

```

```{r}
AIC(var_model)
```

```{r}
forecast_var <- forecast(var_model, h = 18)
plot(forecast_values)
```

```{r}
acf(resid(var_model)[, 1])
acf(resid(var_model)[, 2])
```


###GARCH

```{r}
fit_garch_model <- function(ts_data) {
  garch_fit <- garch(ts_data, order = c(1, 1))

  # Residuals can be considered as errors
  residuals <- garch_fit$residuals
  mae <- mean(abs(residuals), na.rm = TRUE)
  rmse <- sqrt(mean(residuals^2, na.rm = TRUE))

  # AICc calculation
  aicc <- AIC(garch_fit) + (2 * garch_fit$n.coeff * (garch_fit$n.coeff + 1)) / (length(ts_data) - garch_fit$n.coeff - 1)

  return(list(mae = mae, rmse = rmse, aicc = aicc))
}
```


```{r}
fit_garch_model(ts_water)
```

```{r}
ets_model <- ets(ts_water, lambda = "auto")
summary(ets_model)
```



### CROSS VALIDATION



```{r}
k <- 150 # minimum data length for fitting a model
n <- length(ts_vol)

p <- 12 # period/seasonality
H <- 12 # Forecast Horizon

st <- tsp(ts_vol)[1]+(k-2)/p

mae_1 <- matrix(NA,n-k,H)
mae_2 <- matrix(NA,n-k,H)


# Arima expanding window
mae_arima_expand <- matrix(NA,n-k,H)
rmse_arima_expand <- matrix(NA,n-k,H)
aic_arima_expand <- numeric(80)

# Arima sliding window
mae_arima_slide <- matrix(NA,n-k,H)
rmse_arima_slide <- matrix(NA,n-k,H)
aic_arima_slide <- numeric(80)

# ETS expanding window 
mae_ets_expand <- matrix(NA,n-k,H)
rmse_ets_expand <- matrix(NA,n-k,H)
aic_ets_expand <- numeric(80)

# ETS sliding window
mae_ets_slide <- matrix(NA,n-k,H)
rmse_ets_slide <- matrix(NA,n-k,H)
aic_ets_slide <- numeric(80)

for(i in 1:(n-k)) {
    # Expanding Window 
    train_1 <- window(ts_water, end=st + i/p)
    
    # Sliding Window 3,0,1)(1,1,2)[12] 
    train_2 <- window(ts_water, start=st+(i-k+1)/p, end=st+i/p)
    
    test <- window(ts_water, start=st + (i+1)/p, end=st + (i+H)/p)
  
    # Arima - Expanding Window
    fit_1 <- Arima(train_1, order = c(3, 0, 1), seasonal = list(order = c(1, 1, 2), period = 12), lambda = "auto", method='ML')
    fcast_1 <- forecast(fit_1, h=H)
  
    # Arima - Sliding Window
    fit_2 <- Arima(train_2, order = c(3, 0, 1), seasonal = list(order = c(1, 1, 2), period = 12), lambda = "auto", method='ML')
    fcast_2 <- forecast(fit_2, h=H)
  
    # ARIMA errors and AICc - Expanding Window
    mae_arima_expand[i,1:length(test)] <- abs(fcast_1[['mean']]-test)
    rmse_arima_expand[i,1:length(test)] <- (fcast_1[['mean']]-test)^2
    aic_arima_expand[i] <- fit_1$aicc
  
    # ARIMA errors and AICc - Sliding Window
    mae_arima_slide[i,1:length(test)] <- abs(fcast_2[['mean']]-test)
    rmse_arima_slide[i,1:length(test)] <- (fcast_2[['mean']]-test)^2
    aic_arima_slide[i] <- fit_2$aicc
  
    # ARIMAX - Expanding Window
    fit_ets_1 <- ets(train_1, model = "ANA")
    fcast_ets_1 <- forecast(fit_ets_1, h=H)
  
    # ARIMAX - Sliding Window
    fit_ets_2 <- ets(train_2, model = "ANA")
    fcast_ets_2 <- forecast(fit_ets_2, h=H)
  
    # ETS errors and AICc - Expanding Window
    mae_ets_expand[i,1:length(test)] <- abs(fcast_ets_1[['mean']]-test)
    rmse_ets_expand[i,1:length(test)] <- (fcast_ets_1[['mean']]-test)^2
    aic_ets_expand[i] <- fit_ets_1$aicc
  
    # ETS errors and AICc - Sliding Window
    mae_ets_slide[i,1:length(test)] <- abs(fcast_ets_2[['mean']]-test)
    rmse_ets_slide[i,1:length(test)] <- (fcast_ets_2[['mean']]-test)^2
    aic_ets_slide[i] <- fit_ets_2$aicc
}

# Calculate RMSE
rmse_arima_expand <- sqrt(apply(rmse_arima_expand, 2, mean, na.rm = TRUE))
rmse_ets_expand <- sqrt(apply(rmse_ets_expand, 2, mean, na.rm = TRUE))
rmse_arima_slide <- sqrt(apply(rmse_arima_slide, 2, mean, na.rm = TRUE))
rmse_ets_slide <- sqrt(apply(rmse_ets_slide, 2, mean, na.rm = TRUE))

```


```{r}
matplot(cbind(rmse_arima_expand, rmse_arima_slide, rmse_ets_expand, rmse_ets_slide), 
        type = "l", lty = 1, xlab = "Horizon", ylab = "RMSE",
        col = c("red", "blue", "green", "purple"))
legend("topleft", legend = c("ARIMA Expanding", "ARIMA Sliding", "ETS Expanding", "ETS Sliding"), 
       lty = 1, col = c("red", "blue", "green", "purple"), cex = 0.7)
```


```{r}
mae_arima_expand <- colMeans(abs(mae_arima_expand), na.rm=TRUE)
mae_arima_slide <- colMeans(abs(mae_arima_slide), na.rm=TRUE)
mae_ets_expand <- colMeans(abs(mae_ets_expand), na.rm=TRUE)
mae_ets_slide <- colMeans(abs(mae_ets_slide), na.rm=TRUE)

matplot(cbind(mae_arima_expand, mae_arima_slide, mae_ets_expand, mae_ets_slide), 
        type = "l", lty = 1, xlab = "Horizon", ylab = "MAE",
        col = c("red", "blue", "green", "purple"))
legend("topleft", legend = c("ARIMA Expanding", "ARIMA Sliding", "ETS Expanding", "ETS Sliding"),
       lty = 1, col = c("red", "blue", "green", "purple"), cex=0.7)
```


```{r}
plot(aic_arima_expand, type='l', ylab='AICc', xlab='Iteration', main='ARIMA Expanding vs Slide')
lines(aic_arima_slide, col='blue')
legend("topright", legend=c("Expanding", "Sliding"), col=c("black", "blue"), lty = 1)

# ETS
plot(aic_ets_expand, type='l', ylab='AICc', xlab='Iteration', main='ETS Expanding vs Slide')
lines(aic_ets_slide, col='blue')
legend("topright", legend=c("Expanding", "Sliding"), col=c("black", "blue"), lty = 1)
```