---
title: "Temperature-Humidity Index @ Kaunas"
output: 
  flexdashboard::flex_dashboard:
    theme: journal
    orientation: rows
    source_code: embed
    vertical_layout: fill
---

```{r setup, include=FALSE}

# Install / load required R libraries
if (!"pacman" %in% rownames(installed.packages())) install.packages("pacman")
pacman::p_load(xts, prophet, TSrepr, plotly, dygraphs)

# Helper function to estimate various forecasting errors
forecastingErrors <- function(realVals,foreVals,inSamp) {
  RMSE <- TSrepr::rmse(realVals,foreVals)
  MDAE <- TSrepr::mdae(realVals,foreVals)
  MASE <- TSrepr::mae(realVals,foreVals) / inSamp
  MAPE <- TSrepr::mape(realVals,foreVals)
  SMAPE <- TSrepr::smape(realVals,foreVals)
  MAAPE <- TSrepr::maape(realVals,foreVals)
  data.frame(RMSE,MDAE,MASE,MAPE,SMAPE,MAAPE)
}

# Forecasting horizon length
forecastingHorizon <- 28

# Load prepared xts data
myDataFull <- readRDS('waqi-covid19-airqualitydata.rds')
myDataTrain <- myDataFull[1:(nrow(myDataFull)-forecastingHorizon-1)]

```

```{r backtesting, include=FALSE}

# Denominator for MASE forecasting error calculations
inSampleMAE <- mean(abs(diff(myDataTrain$THI)), na.rm=T)

# Univariate forecasting of temperature for myDataTrain
df <- data.frame(ds=as.Date(time(myDataTrain)), y=as.numeric(myDataTrain[,'temperature']))
tempModel <- prophet(df, yearly.seasonality=T, weekly.seasonality=T, daily.seasonality=F)
myFuture <- make_future_dataframe(tempModel, periods = forecastingHorizon, freq='day')
tempForecast <- predict(tempModel, myFuture)

# Univariate forecasting of humidity for myDataTrain
df <- data.frame(ds=as.Date(time(myDataTrain)), y=as.numeric(myDataTrain[,'humidity']))
humidModel <- prophet(df, yearly.seasonality=T, weekly.seasonality=T, daily.seasonality=F)
myFuture <- make_future_dataframe(humidModel, periods = forecastingHorizon, freq='day')
humidForecast <- predict(humidModel, myFuture)

# Multivariate forecasting of THI for myDataTrain
df <- data.frame(ds=as.Date(time(myDataTrain)), y=as.numeric(myDataTrain[,'THI']),
                 myDataTrain[,'temperature'], myDataTrain[,'humidity'])
indexModel <- prophet(yearly.seasonality=T, weekly.seasonality=T, daily.seasonality=F)
indexModel <- add_regressor(indexModel, 'temperature')
indexModel <- add_regressor(indexModel, 'humidity')
indexModel <- fit.prophet(indexModel, df)
myFuture <- make_future_dataframe(indexModel, periods = forecastingHorizon, freq='day')
myFuture <- cbind(myFuture, temperature=tempForecast$yhat, humidity=humidForecast$yhat)
indexForecast <- predict(indexModel, myFuture)

# Calculate forecasting errors for THI on testing data (last observations)
tempPred <- tail(tempForecast$yhat, forecastingHorizon)
humidPred <- tail(humidForecast$yhat, forecastingHorizon)
THI_uniformula <- c(rep(NA, 2*forecastingHorizon), 0.8*tempPred + (humidPred/100)*(tempPred-14.4) + 46.4)
THI_multiforec <- c(rep(NA, 2*forecastingHorizon), tail(indexForecast$yhat, forecastingHorizon))
forecastedTHI <- cbind(tail(myDataFull$THI, 3*forecastingHorizon), THI_uniformula, THI_multiforec)
idx <- 1:forecastingHorizon + 2*forecastingHorizon
errorsTable <- rbind(THI_uniformula=forecastingErrors(forecastedTHI[idx,1], forecastedTHI[idx,2], inSampleMAE),
                     THI_multiforec=forecastingErrors(forecastedTHI[idx,1], forecastedTHI[idx,3], inSampleMAE))

```

```{r forecasting, include=FALSE}

# Univariate forecasting of temperature for myDataFull
df <- data.frame(ds=as.Date(time(myDataFull)), y=as.numeric(myDataFull[,'temperature']))
tempModel <- prophet(df, yearly.seasonality=T, weekly.seasonality=T, daily.seasonality=F)
myFuture <- make_future_dataframe(tempModel, periods = forecastingHorizon, freq='day')
tempForecast <- predict(tempModel, myFuture)

# Univariate forecasting of humidity for myDataFull
df <- data.frame(ds=as.Date(time(myDataFull)), y=as.numeric(myDataFull[,'humidity']))
humidModel <- prophet(df, yearly.seasonality=T, weekly.seasonality=T, daily.seasonality=F)
myFuture <- make_future_dataframe(humidModel, periods = forecastingHorizon, freq='day')
humidForecast <- predict(humidModel, myFuture)

# Multivariate forecasting of THI for myDataFull
df <- data.frame(ds=as.Date(time(myDataFull)), y=as.numeric(myDataFull[,'THI']),
                 myDataFull[,'temperature'], myDataFull[,'humidity'])
indexModel <- prophet(yearly.seasonality=T, weekly.seasonality=T, daily.seasonality=F)
indexModel <- add_regressor(indexModel, 'temperature')
indexModel <- add_regressor(indexModel, 'humidity')
indexModel <- fit.prophet(indexModel, df)
myFuture <- make_future_dataframe(indexModel, periods = forecastingHorizon, freq='day')
myFuture <- cbind(myFuture, temperature=tempForecast$yhat, humidity=humidForecast$yhat)
indexForecast <- predict(indexModel, myFuture)

```

Backtesting
=======================================================================

Row
-------------------------------------

### **Reality (black) vs. forecast**

```{r plot.xts, fig.width = 10}

gg <- xts::plot.xts(forecastedTHI)
xts::addLegend(legend.names = names(forecastedTHI), lty=1, col=1:ncol(forecastedTHI))

```

Row {data-height=110}
-------------------------------------

### **Table of forecasting errors**

```{r table}

knitr::kable(errorsTable, digits=3)

```

Components
=======================================================================

```{r prophet, fig.show='hide'}

gg <- prophet_plot_components(indexModel, indexForecast)

```

Row {.tabset .tabset-fade}
-------------------------------------

### **Trend component**

```{r prophet1}

ggplotly(gg[[1]])

```

### **Extra regressors**

```{r prophet4}

ggplotly(gg[[4]])

```

Row {.tabset .tabset-fade}
-------------------------------------

### **Yearly seasonality**

```{r prophet3}

ggplotly(gg[[3]])

```

### **Weekly seasonality**

```{r prophet2}

ggplotly(gg[[2]])

```

Forecasts (ggplotly)
=======================================================================

Row
-------------------------------------

### **Temperature**

```{r ggplotly-temp}

ggplotly(plot(tempModel, tempForecast, ylabel='Celsius (°C)', xlabel='Time (daily)'))

```

Row
-------------------------------------

### **Humidity**

```{r ggplotly-humid}

ggplotly(plot(humidModel, humidForecast, ylabel='Percentage (%)', xlabel='Time (daily)'))

```

Row
-------------------------------------

### **Temperature-Humidity Index**

```{r ggplotly-THI}

ggplotly(plot(indexModel, indexForecast, ylabel='THI', xlabel='Time (daily)'))

```

Forecasts (dygraphs)
=======================================================================

```{r dygraphs, include=FALSE}

timeStamps <- tail(time(myDataFull),1) + 1:forecastingHorizon

```

Row
-------------------------------------

### **Temperature (°C)**

```{r dygraphs-temp}

df <- tempForecast[tempForecast$ds %in% timeStamps, c('yhat_lower','yhat','yhat_upper')]
names(df) <- c('Lo.80', 'Point.Forecast', 'Hi.80')
myFore <- xts(df, order.by = timeStamps)
plotData <- cbind(myDataFull[,'temperature'], myFore)
dygraph(plotData, group = "Air-Data") %>%
  dySeries(c('temperature'), label = 'Temperature') %>%
  dySeries(c("Lo.80", "Point.Forecast", "Hi.80"), label = "Forecast") %>%
  dyCrosshair(direction = "vertical")

```

Row
-------------------------------------

### **Humidity (%)**

```{r dygraphs-humid}

df <- humidForecast[humidForecast$ds %in% timeStamps, c('yhat_lower','yhat','yhat_upper')]
names(df) <- c('Lo.80', 'Point.Forecast', 'Hi.80')
myFore <- xts(df, order.by = timeStamps)
plotData <- cbind(myDataFull[,'humidity'], myFore)
dygraph(plotData, group = "Air-Data") %>%
  dySeries(c('humidity'), label = 'Humidity') %>%
  dySeries(c("Lo.80", "Point.Forecast", "Hi.80"), label = "Forecast") %>%
  dyCrosshair(direction = "vertical")

```

Row
-------------------------------------

### **Temperature-Humidity Index**

```{r dygraphs-THI}

df <- indexForecast[indexForecast$ds %in% timeStamps, c('yhat_lower','yhat','yhat_upper')]
names(df) <- c('Lo.80', 'Point.Forecast', 'Hi.80')
myFore <- xts(df, order.by = timeStamps)
plotData <- cbind(myDataFull[,'THI'], myFore)
dygraph(plotData, group = "Air-Data") %>%
  dySeries(c('THI'), label = 'THI') %>%
  dySeries(c("Lo.80", "Point.Forecast", "Hi.80"), label = "Forecast") %>%
  dyRangeSelector() %>% dyCrosshair(direction = "vertical")

```