---
title: "PosouzenÃ­ moÅ¾nostÃ­ vyuÅ¾itÃ­ rÅ¯znÃ½ch predikÄnÃ­ch modelÅ¯ pro dostupnÃ¡ dopravnÃ­ data"
---
Tento Quarto markdown dokument obsahuje kÃ³d, pomocÃ­ nÄhoÅ¾ byla vypracovÃ¡na ExplorativnÃ­ datovÃ¡ analÃ½za, tvorba predikÄnÃ­ch modelÅ¯ a validace vytvoÅenÃ½ch modelÅ¯. Pouze ÄÃ¡st kÃ³du byla pouÅ¾ita v pÅedklÃ¡danÃ© prÃ¡ci.

```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(janitor)
library(magrittr)
library(ggsci)
library(lubridate)
library(leaflet)
library(XML)
library(xml2)
library(zoo)

library(xgboost)
library(broom)
library(tidymodels)
library(modeltime)
library(modeltime.resample)
library(timetk)
library(tsfknn)

Sys.setlocale(category = "LC_TIME", locale = "cs_CZ.UTF-8")
```

```{r message=FALSE, warning=FALSE}
senzor_data_colnames <- c("IdDetektor", "DatumCas", "Intenzita", "IntenzitaN", "Obsazenost", "Rychlost", "Stav", "TypVozidla", "Trvani100", "RychlostHistorie", "TypVozidla10")

#nacteni dat z stacionarniho detektoru
senzor_data <- read_delim("DOPR_D_20210126.csv", delim = "|", col_names = senzor_data_colnames) %>%
  janitor::clean_names()
  
```

# Explorativni Datova Analyza

## Dobrichovice

```{r}
dobrichovice_coords <- tibble(lat = c(49.929269, 49.953629), long = c(14.273824, 14.319726))

leaflet() %>%
  addTiles() %>%
  addPolylines(data = dobrichovice_coords, lng = ~long, lat = ~lat)
```

## Stacionarni detektor

| NÃ¡zev sloupce     | Popis sloupce                                                         |
|-------------------|-----------------------------------------------------|
| id_detektor       | identifikaÄnÃ­ ÄÃ­slo detektoru s rozliÅ¡enÃ­m smÄru a pruhu              |
| datum_cas         |                                                                       |
| intenzita         | poÄet vozidel ze ÄasovÃ½ interval                                      |
| intenzita_n       | normovanÃ¡ intenzita - motorka mÃ¡ menÅ¡Ã­ vÃ¡hu neÅ¾ nÃ¡klaÄÃ¡k              |
| obsazenost        | obszenost detektoru v procentech                                      |
| rychlost          | rychlost vozidla v kmh                                                |
| stav              | 1- funkÄnÃ­, \<0 nefunkÄnÃ­ (moc nedava smysl)                          |
| typ_vozidla       | zjednoduÅ¡enÃ½ typ vozidla, 1 - bike, 2 - car, 3 - van, 4 - truck       |
| trvani100         | ÄasovÃ½ rozestup mezi daty v historii rychlosti ve stovkÃ¡ch milisekund |
| rychlost_historie | historie rychlostÃ­ k danÃ© detekci vozidla                             |
| typ_vozidla10     | pÅesnÄjÅ¡Ã­ detekce typu vozidla                                        |

```{r}

senzor_data %<>%
  separate(datum_cas, into = c("datum", "cas"), sep = " ") %>%
  mutate(datum = as_date(datum)) %>%
  #mutate(cas = hms(cas)) %>%
  mutate(typ_vozidla = factor(typ_vozidla)) %>%
  mutate(typ_vozidla = recode_factor(typ_vozidla,"1" = "bike", "2" =  "car", "3" = "van", "4" = "truck")) %>%
  mutate(id_detektor = as_factor(id_detektor)) %>%
  mutate(
    detektor_lokace = case_when(id_detektor == "42057111" ~ "z DobÅichovic do ÄernoÅ¡ic",
                                id_detektor == "42057112" ~ "z ÄernoÅ¡ic do DobÅichovic", 
                                id_detektor == "42057211" ~ "z DobÅichovic do Let",
                                id_detektor == "42057212" ~ "z Let do DobÅichovic" )
    )

```

```{r}
senzor_data %>%
    summary()
```

### Obecne info o celem datovem souboru

-   2 136 407 zaznamu
-   Nejstarsi datum 2020-06-19
-   Nejnovejsi datum 2021-01-26

### Lokace detektoru podle ID

-   42057111 - ulice PraÅ¾skÃ¡, z DobÅichovic do ÄernoÅ¡ic
-   42057112 - ulice PraÅ¾skÃ¡, z ÄernoÅ¡ic do DobÅichovic
-   42057211 - ulice PraÅ¾skÃ¡, z DobÅichovic do Let
-   42057212 - ulice PraÅ¾skÃ¡, z Let do DobÅichovic

```{r}
senzor_data %>%
  head()
```

```{r}
senzor_data %>%
  group_by(typ_vozidla) %>%
  summarise(count = n()) %>%
  mutate(freq = round(count / sum(count), 2))
```

### VÅ¡echna data ze vÅ¡ech senzorÅ¯

```{r}
senzor_data %>%
ggplot(aes(x = rychlost)) +
  geom_histogram(binwidth = 1, fill = "#B24745FF") +
  theme_classic() +
  labs(title = "Histogram rychlosti vozidel ze vÅ¡ech detektorÅ¯", x = "Rychlost [km/h]" , y = "PoÄet vozidel [-]")

ggsave("str_dobr_all_hist_speed.svg")
```

-   Hodne outlieru - dlouha osa x

-   Hodne mereni rychlosti 0

    -   U mereni, kde rychlost = 0 vzdy chybi sloupec rychlost_historie

    -   Rychlost = 0 u 350 174 zaznamu

```{r}
senzor_data %>%
  filter(rychlost > 0 & rychlost < 70) %>%
ggplot(aes(x = rychlost)) +
  geom_histogram(binwidth = 1, fill = "#B24745FF") +
  scale_fill_jama() +
  theme_classic() +
  facet_wrap(vars(detektor_lokace)) +
  scale_x_continuous(limits = c(0, 70), breaks = seq(0, 70, 10)) +
  labs(title = "Histogram rychlosti vozidel ze vÅ¡ech detektorÅ¯", subtitle = "RozdÄleno podle lokece detektoru", x = "Rychlost [km/h]" , y = "PoÄet vozidel [-]") 

ggsave("str_dobr_all_hist_speed_facet.svg")
```

```{r}
senzor_data %>%
    filter(rychlost > 0 & rychlost < 70) %>%
ggplot(aes(x = rychlost)) +
  geom_histogram(binwidth = 1, fill = "#B24745FF") +
  theme_classic() +
  labs(title = "Histogram rychlosti vozidel ze vÅ¡ech senzorÅ¯",
       subtitle = "Pouze rychlost menÅ¡Ã­ neÅ¾ 70 a bez mÄÅenÃ­ kde rychlost rovna nule",
       x = "Rychlost [km/h]" ,
       y = "PoÄet vozidel [-]")
```

```{r}
senzor_data %>%
    filter(rychlost > 70) %>%
ggplot(aes(x = rychlost)) +
  geom_histogram(binwidth = 1, fill = "#B24745FF") +
  theme_classic() +
  labs(title = "Histogram rychlosti vozidel ze vÅ¡ech senzorÅ¯",
       subtitle = "Pouze mÄÅenÃ­ kde rychlost pÅesÃ¡hla 70 km/h",
       x = "Rychlost [km/h]" ,
       y = "PoÄet vozidel [-]")
```

```{r}
senzor_data %>%
  #filter(datum == "2020-12-01") %>%
  filter(typ_vozidla == "car") %>%
  mutate(fake_datum = ymd_hms(paste("2020-01-05",cas))) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  relocate(datum, .after = cas) %>% 
  mutate(hodina = hour(datum)) %>%
  #group_by(hour(datum), minute(datum)) %>%
  group_by(cas_mereni =  round_date(fake_datum, "5 minutes")) %>%
    summarise(n_vozidel = n(), prumerna_rychlost = mean(rychlost), sd_rychlost = sd(rychlost)) %>%
    #separate(cas_mereni, into = c("datum", "cas"), sep = " ") %>%
    #select(-datum) %>%
    #mutate(cas = hms(cas)) %>%
ggplot(aes(x = cas_mereni, y = n_vozidel)) +
  geom_line() +
  scale_x_datetime(date_labels = "%H:%M") +
  theme_classic() +
  labs(title = "Celkovy pocet zmerenych vozidel kazdych 5 minut na vsech detektorech")
```

```{r}
senzor_data %>%
  #filter(datum == "2020-12-01") %>%
  filter(typ_vozidla == "car") %>%
  mutate(fake_datum = ymd_hms(paste("2020-01-05",cas))) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  relocate(datum, .after = cas) %>% 
  mutate(hodina = hour(datum)) %>%
  #group_by(hour(datum), minute(datum)) %>%
  group_by(detektor_lokace, den_mereni = date(datum), cas_mereni =  round_date(fake_datum, "5 minutes")) %>%
    summarise(n_vozidel = n(), prumerna_rychlost = mean(rychlost), sd_rychlost = sd(rychlost)) %>%
    mutate(intenzita = n_vozidel/(5/60)) %>%
  ungroup() %>%
ggplot(aes(x=cas_mereni, y = intenzita, group = den_mereni)) +
  geom_line(alpha = 0.01) +
  facet_wrap(vars(detektor_lokace)) +
  theme_classic() +
  scale_x_datetime(date_labels = "%H:%M") +
  labs(x = "Äas", y = "Intenzita [voz/hod]", title = "DennÃ­ dopravnÃ­ intenzita", subtitle = "RozdÄleno podle lokece detektoru")

ggsave("str_dobr_all_intensity.svg")
```

```{r}
senzor_data %>%
  #filter(datum == "2020-12-01") %>%
  filter(typ_vozidla == "car") %>%
  filter(rychlost != 0) %>%
  mutate(fake_datum = ymd_hms(paste("2020-01-05",cas))) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  relocate(datum, .after = cas) %>% 
  mutate(hodina = hour(datum)) %>%
  #group_by(hour(datum), minute(datum)) %>%
  group_by(detektor_lokace, den_mereni = date(datum), cas_mereni =  round_date(fake_datum, "5 minutes")) %>%
    summarise(n_vozidel = n(), prumerna_rychlost = mean(rychlost), sd_rychlost = sd(rychlost)) %>%
ggplot(aes(x=cas_mereni, y = prumerna_rychlost, group = den_mereni, color = "darkred")) +
  geom_line(alpha = 0.01) +
  facet_wrap(vars(detektor_lokace)) +
  theme_classic() +
  scale_x_datetime(date_labels = "%H:%M") +
  theme(legend.position = "none") +
  labs(x = "Äas", y = "Rychlost [km/h]", title = "DennÃ­ rychlosti vozidel", subtitle = "RozdÄleno podle lokece detektoru")

#ggsave("str_dobr_all_intensity.svg")
```

```{r}
senzor_data %>%
  #filter(datum == "2020-12-01") %>%
  filter(typ_vozidla == "car") %>%
  mutate(fake_datum = ymd_hms(paste("2020-01-05",cas))) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  relocate(datum, .after = cas) %>% 
  mutate(hodina = hour(datum)) %>%
  #group_by(hour(datum), minute(datum)) %>%
  group_by(detektor_lokace, den_mereni = date(datum)) %>%
    summarise(n_vozidel = n(), prumerna_rychlost = mean(rychlost), sd_rychlost = sd(rychlost)) %>%
ggplot(aes(x=den_mereni, y = n_vozidel)) +
  geom_step() +
  facet_wrap(vars(detektor_lokace)) +
  theme_classic() +
  labs(title = "DennÃ­ intenzity", x = "Datum", y = "Intenzita [voz/den]")

ggsave("str_dobr_daily_intensity.svg")
```

### Dobrichovice - Cernosice

-   Celkem 500 668 zaznamu (vozidel)

```{r}
senzor_data %>%
  filter(id_detektor == 42057111) %>%
  filter(typ_vozidla == "car") %>%
  group_by(den_mereni =  round_date(datum, "1 day")) %>%
    summarise(n_vozidel = n(), prumerna_rychlost = mean(rychlost), sd_rychlost = sd(rychlost)) %>%
  ungroup() %>%
ggplot(aes(x = den_mereni, y = n_vozidel)) +
  geom_step() +
  #geom_bar(stat = "identity") +
  theme_classic() +
  labs(title = "PoÄet namÄÅenÃ½ch vozidel",
       subtitle = "Na detektoru ve smÄru z Dobrichovic do Cernosic", 
       x = "Datum",
       y = "PoÄet vozidel")
```

```{r}
senzor_data %>%
  #filter(id_detektor == 42057111) %>%
  filter(typ_vozidla == "car") %>%
  filter(datum == "2020-09-20") %>%
  mutate(fake_datum = ymd_hms(paste("2020-09-20",cas))) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  relocate(datum, .after = cas) %>% 
  mutate(hodina = hour(datum)) %>%
  #group_by(hour(datum), minute(datum)) %>%
  group_by(detektor_lokace, cas_mereni =  round_date(fake_datum, "15 minutes")) %>%
    summarise(n_vozidel = n(), prumerna_rychlost = mean(rychlost), sd_rychlost = sd(rychlost)) %>%
    mutate(intenzita = n_vozidel/(15/60)) %>%
  ggplot(aes(x = cas_mereni, y = intenzita)) +
  geom_line() +
  scale_x_datetime(date_labels = "%H:%M") +
  theme_classic() +
  facet_wrap(vars(detektor_lokace)) +
  labs(title = "DopravnÃ­ intenzita dne 20. 9. 2020",
       subtitle = "RozdÄleno dle lokace",
       x = "Äas",
       y = "Intenzita [voz/hod]")

ggsave("str_dobr_day_intensity.svg")
```

```{r}
senzor_data %>%
  #filter(id_detektor == 42057111) %>%
  filter(rychlost != 0) %>%
  filter(typ_vozidla == "car") %>%
  filter(datum == "2020-09-20") %>%
  mutate(fake_datum = ymd_hms(paste("2020-09-20",cas))) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  relocate(datum, .after = cas) %>% 
  mutate(hodina = hour(datum)) %>%
  #group_by(hour(datum), minute(datum)) %>%
  group_by(detektor_lokace, cas_mereni =  round_date(fake_datum, "15 minutes")) %>%
    summarise(n_vozidel = n(), prumerna_rychlost = mean(rychlost), sd_rychlost = sd(rychlost)) %>%
  ggplot(aes(x = cas_mereni, y = prumerna_rychlost)) +
  geom_line(color = "#B24745FF") +
  geom_ribbon(aes(y = prumerna_rychlost, ymin = prumerna_rychlost - sd_rychlost, ymax = prumerna_rychlost + sd_rychlost),fill = "red", alpha = .2) +
  scale_x_datetime(date_labels = "%H:%M") +
  theme_classic() +
  facet_wrap(vars(detektor_lokace)) +
  labs(title = "Rychlost vozidel",
       subtitle = "Dne 20. 9. 2020, rozdÄleno dle lokace",
       x = "Äas",
       y = "Rychlost [km/h]")

ggsave("str_dobr_day_speed.svg")
```

```{r}
senzor_data %>%
  filter(datum == "2020-12-01") %>%
  filter(typ_vozidla == "car") %>%
  mutate(whole_date = ymd_hms(paste(datum,cas))) %>%
  relocate(whole_date, .after = cas)
```

## FCD

```{r}
fcd_data <- read_csv2("11_2019.csv") %>%
  janitor::clean_names()
```

```{r}
fcd_data %<>%
  filter(smer %in% c('z DobÅichovice smÄr ÄernoÅ¡ice', 'z ÄernoÅ¡ice smÄr DobÅichovice', 'z Lety smÄr DobÅichovice', 'z DobÅichovice smÄr Lety')) %>%
  mutate(smer = str_replace(smer, "smÄr", "do")) %>%
  mutate(smer = str_replace(smer, "ÄernoÅ¡ice", "ÄernoÅ¡ic")) %>%
  mutate(smer = str_replace(smer, "DobÅichovice", "DobÅichovic")) %>%
  mutate(smer = str_replace(smer, "Lety", "Let")) %>%
  mutate(datum = str_remove(datum, "\\s\\+\\d{2}:\\d{2}")) %>%
  separate(datum, into = c("datum", "cas"), sep = " ") %>%
  mutate(datum = dmy(datum)) %>%
  rename(n_osobni_vozidla = vd1_1_aktualni_pocet_plovoucich_vozidel_osobnich) %>%
  rename(n_nakladni_vozidla = vd1_2_aktualni_pocet_plovoucich_vozidel_nakladnich) %>%
  mutate(pocet_vozidel = n_osobni_vozidla + n_nakladni_vozidla) %>%
  rename(rychlost = vd2_aktualni_vypoctena_rychlost_dopravniho_proudu_km_h) %>%
  rename(cas_prujezdu = vd3_aktualni_cas_prujezdu_definovaneho_tmc_segmentu_s) %>%
  rename(zpozdeni = vd4_aktualni_zpozdeni_na_definovanem_segmentu_s) %>%
  rename(lokalizace_kolony = vd7_1_lokalizace_kolony_na_segmentu) %>%
  rename(delka_kolony = vd7_2_delka_kolony_m) %>%
  rename(vzdalenost_kolony_od_pocatku_tmc = vd7_3_vzdalenost_cela_kolony_od_pocatku_segmentu_m) %>%
  rename(stupen_dopravy = vd10_parametr_level_of_service_lo_s) #%>%
  #select(-c(oblast, komunikace, smer, tmc_id, vd5_typicka_prumerna_rychlost_pri_volnem_prujezdu_segmentu_km_h, vd6_doba_volneho_prujezdu_s, vd8_mira_spolehlivosti, vd9_reakcni_doby_systemu_s))
```

```{r}
unique(fcd_data$smer)
```

```{r}
fcd_data %>%
  group_by(smer, den_mereni =  round_date(datum, "1 day")) %>%
    summarise(n_vozidel = sum(pocet_vozidel), prumerna_rychlost = mean(rychlost), sd_rychlost = sd(rychlost)) %>%
ggplot(aes(x = den_mereni, y = n_vozidel)) +
  geom_step() +
  #geom_bar(stat = "identity") +
  theme_classic() +
  scale_x_date(date_labels = "%d. %m") +
  facet_wrap(vars(smer)) +
  labs(title = "DennÃ­ poÄty vozidel", subtitle = "Z FCD", y = "PoÄet vozidel [-]", x = "Datum")

ggsave("fcd_dobr_daily_intensity.svg")
```

```{r}
fcd_data %>%
  filter(datum == '2019-11-15') %>%
  mutate(datum = as.character(datum)) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  select(-cas) %>%
ggplot(aes(x = datum, y = pocet_vozidel)) +
  geom_line() +
  scale_x_datetime(date_labels = "%H:%M") +
  theme_classic() +
  labs(title = "PoÄet vozidel", x = "Äas", y = "PoÄet vozidel")
```

```{r}
fcd_data %>%
  filter(datum == '2019-11-15') %>%
  mutate(datum = as.character(datum)) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  select(-cas) %>%
  group_by(smer, cas_mereni =  round_date(datum, "15 minutes")) %>%
    summarise(n_vozidel = n(), prumerna_rychlost = mean(rychlost), sd_rychlost = sd(rychlost)) %>%
  ungroup() %>%
ggplot(aes(x = cas_mereni, y = n_vozidel)) +
  geom_line() +
  scale_x_datetime(date_labels = "%H:%M") +
  scale_y_continuous(limits = c(0,16)) +
  theme_classic() +
  facet_wrap(vars(smer)) +
  labs(title = "PoÄet vozidel z FCD", subtitle = "V jeden nÃ¡hodnÃ½ den", x = "Äas", y = "PoÄet vozidel [-]")

ggsave("fcd_dobr_day_intensity.svg")
```

```{r}
fcd_data %>%
  filter(datum == '2019-11-15') %>%
  mutate(datum = as.character(datum)) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  select(-cas) %>%
  group_by(smer, cas_mereni = round_date(datum, "10 minutes")) %>%
    summarise(n_vozidel = n(), prumerna_rychlost = mean(rychlost), sd_rychlost = sd(rychlost)) %>%
  ungroup() %>%
ggplot(aes(x = cas_mereni, y = prumerna_rychlost)) +
  geom_line(color = "red") +
  geom_ribbon(aes(y = prumerna_rychlost, ymin = prumerna_rychlost - sd_rychlost, ymax = prumerna_rychlost + sd_rychlost),fill = "red", alpha = .2) +
  scale_x_datetime(date_labels = "%H:%M") +
  theme_classic() +
  facet_wrap(vars(smer)) +
  labs(title = "Rychlost vozidel z FCD", x = "Äas", y = "Rychlost [km/h]")

ggsave("fcd_dobr_day_speed.svg")
```

```{r}
fcd_data %>%
  filter(datum == '2019-11-15') %>%
  mutate(datum = as.character(datum)) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  select(-cas) %>%
ggplot(aes(x = datum, y = zpozdeni)) +
  geom_line() +
  theme_classic() +
  scale_x_datetime(date_labels = "%H:%M") +
  labs(title = "ZpoÅ¾dÄnÃ­", x = "Äas", y = "ZpoÅ¾dÄnÃ­ [s]")
```

```{r}
fcd_data %>%
  #filter(datum == '2019-11-15') %>%
  mutate(datum = as.character(datum)) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  select(-cas) %>%
ggplot(aes(x = zpozdeni)) +
  geom_histogram(binwidth = 10, fill = "#6A6599FF") + 
  theme_classic() +
  facet_wrap(vars(smer)) +
  labs(title = "Histogram zpoÅ¾dÄnÃ­", x = "ZpoÅ¾dÄnÃ­ [s]", y = "PoÄet")

ggsave("fcd_dobr_daily_delay.svg")
```



```{r}
#nestacionarni data
set.seed(8)
data.frame(
  day = as.Date("2022-01-01") - 0:364,
  value = cumsum(rnorm(n=365, mean = 1, sd = 20))
) %>%
  ggplot(aes(x = day, y = value)) +
    geom_line() +
    theme_classic() +
    scale_x_date(date_breaks = "3 months") +
    theme(text = element_text(size = 20)) +
    labs(x = "Datum", y = "PromÄnnÃ¡")

ggsave("other_nestacionarni_data.svg")

set.seed(8)
data.frame(
  day = as.Date("2022-01-01") - 0:364,
  value = rnorm(n=365, mean = 1, sd = 2.5)
) %>%
  ggplot(aes(x = day, y = value)) +
  geom_line() +
  theme_classic() +
  scale_x_date(date_breaks = "3 months") +
  scale_y_continuous(limits = c(-10,10)) +
  theme(text = element_text(size = 20)) +
  labs(x = "Datum", y = "PromÄnnÃ¡")

ggsave("other_stacionarni_data.svg")
```


## Tvorba modelu

### Data ze senzoru - Dobrichovice

```{r}
senzor_data %>%
  filter(id_detektor == 42057111) %>%
  filter(typ_vozidla == "car") %>%
  filter(cas >= '09:00:00' & cas <= '17:00:00') %>%
  filter(datum == "2020-09-01") %>%
  mutate(fake_datum = ymd_hms(paste("2020-09-01",cas))) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  relocate(datum, .after = cas) %>% 
  mutate(hodina = hour(datum)) %>%
  #group_by(hour(datum), minute(datum)) %>%
  group_by(cas_mereni =  round_date(fake_datum, "30 minutes")) %>%
    summarise(n_vozidel = n(), prumerna_rychlost = mean(rychlost)) %>%
  ungroup()
```

```{r}
d_c_data <- senzor_data %>%
  filter(id_detektor == 42057111) %>%
  filter(typ_vozidla == "car") %>%
  filter(datum >= "2020-09-01" & datum <= "2020-12-04") %>%
  mutate(fake_datum = ymd_hms(paste("2020-09-01",cas))) %>%
  mutate(datum = ymd_hms(paste(datum,cas))) %>%
  relocate(datum, .after = cas) %>% 
  mutate(hodina = hour(datum)) %>%
  group_by(cas_mereni =  round_date(datum, "5 minutes")) %>%
    summarise(n_vozidel = n(), prumerna_rychlost = mean(rychlost)) %>%
    mutate(intenzita = n_vozidel/(5/60)) %>%
  ungroup() %>%
  select(-n_vozidel) %>%
  mutate(cas_mereni = as_datetime(cas_mereni)) #%>%
  #select(-prumerna_rychlost)
  

d_c_data %>%
  plot_time_series(cas_mereni, intenzita, .interactive = TRUE, .smooth = FALSE)
```

```{r}
d_c_data_splits <- d_c_data %>%
  time_series_split(date_var = cas_mereni, assess = "16 hours", cumulative = TRUE)

d_c_data_splits %>%
  tk_time_series_cv_plan() %>%
  plot_time_series_cv_plan(cas_mereni, intenzita, .interactive = TRUE)
```

```{r}
cv_splits <- d_c_data %>%
  time_series_cv(date_var = cas_mereni, assess = "6 day", initial = "60 days", skip = "5 days", slice_limit = 6)

cv_splits %>% 
  tk_time_series_cv_plan() %>%
  plot_time_series_cv_plan(cas_mereni, intenzita, .interactive = FALSE, .facet_ncol = 2)

model_fit_prophet  <- prophet_reg() %>%
   set_engine('prophet') %>% 
  fit(intenzita ~ cas_mereni, data = training(cv_splits$splits[[1]]))

prophet_tbl <- modeltime_table(model_fit_prophet)

calibration_tbl <- prophet_tbl %>%
  modeltime_calibrate(testing(cv_splits$splits[[1]]))

forecast_panel_tbl <- calibration_tbl %>%
  modeltime_forecast(
    new_data      = testing(cv_splits$splits[[1]]),
    actual_data   = d_c_data,
    keep_data = TRUE
  ) 

resample_results <- prophet_tbl %>%
  modeltime_fit_resamples(
    resamples = cv_splits,
    control   = control_resamples(verbose = FALSE)
  )

resample_results %>%
  modeltime_resample_accuracy(summary_fns = mean) %>%
    table_modeltime_accuracy(.interactive = FALSE)

#-----
# !! FOR LOOP HERE :(
calibration_tbl_test_1 <- prophet_tbl %>%
  modeltime_calibrate(training(cv_splits$splits[[1]]))

calibration_tbl_test_1 %>%
    modeltime_accuracy() %>%
    table_modeltime_accuracy(
        .interactive = FALSE
    )

calibration_tbl_test_2 <- prophet_tbl %>%
  modeltime_calibrate(training(cv_splits$splits[[2]]))

calibration_tbl_test_2 %>%
    modeltime_accuracy() %>%
    table_modeltime_accuracy(
        .interactive = FALSE
    )


prophet_tbl %>%
  modeltime_calibrate(training(cv_splits$splits[[3]])) %>%
    modeltime_accuracy() %>%
    table_modeltime_accuracy(
        .interactive = FALSE
    )
```

```{r}
d_c_data %>%
  time_series_cv(date_var = cas_mereni, assess = "1 day", initial = "4 days", skip = "1 days", slice_limit = 6) %>%
    tk_time_series_cv_plan() %>%
  mutate(.id = str_replace(.id, "Slice", "Vzorek ")) %>%
ggplot(aes(x = cas_mereni, y = intenzita, color = .key)) +
  geom_line() +
  facet_wrap(vars(.id)) +
  scale_color_aaas(labels = c("TestovacÃ­", "TrÃ©novacÃ­")) +
  scale_x_datetime(date_labels = "%d.%m") + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position="bottom") +
  labs(color = "Druh dat", x = "Datum", y = "Intenzita [voz/hod]")

ggsave("cross_validation_example.svg")
```

```{r}

recipe_obj <- recipe(intenzita ~ ., training(d_c_data_splits)) %>%
  step_timeseries_signature(cas_mereni) %>%
  #step_rm(cas_mereni) %>%
  step_zv(all_predictors()) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE)

model_fit_prophet_boost  <- prophet_boost() %>%
    set_engine('prophet_xgboost') %>% 
    fit(intenzita ~ cas_mereni, data = training(d_c_data_splits))

model_fit_prophet  <- prophet_reg() %>%
   set_engine('prophet') %>% 
  fit(intenzita ~ cas_mereni, data = training(d_c_data_splits))

workflow_prophet_boost <- workflow() %>%
  add_model(
     prophet_boost() %>% set_engine('prophet_xgboost') 
  ) %>%
  add_recipe(recipe_obj) %>%
  fit(training(d_c_data_splits))

train_preds <- model_fit_prophet %>% augment(training(d_c_data_splits))

prophet_tbl <- modeltime_table(model_fit_prophet, model_fit_prophet_boost, workflow_prophet_boost)

calibration_tbl <- prophet_tbl %>% 
  modeltime_calibrate(new_data = testing(d_c_data_splits))

forecast_tbl <- calibration_tbl %>%
    modeltime_forecast(
        new_data    = testing(d_c_data_splits),
        actual_data = d_c_data,
        keep_data = TRUE
    )


forecast_tbl %>% plot_modeltime_forecast()

calibration_tbl %>%
    modeltime_accuracy() %>%
    table_modeltime_accuracy(
        .interactive = FALSE
    )


calibration_tbl_test <- prophet_tbl %>% 
  modeltime_calibrate(new_data = training(d_c_data_splits))


forecast_tbl_test <- calibration_tbl_test %>%
    modeltime_forecast(
        new_data    = training(d_c_data_splits),
        actual_data = d_c_data,
        keep_data = TRUE
    )


forecast_tbl_test %>% plot_modeltime_forecast()

calibration_tbl_test %>%
    modeltime_accuracy() %>%
    table_modeltime_accuracy(
        .interactive = FALSE
    )
```

## Prediction models training

### no CV
```{r}
d_c_data_splits <- d_c_data %>%
  select(-prumerna_rychlost) %>%
  time_series_split(date_var = cas_mereni, assess = "16 hours", cumulative = TRUE)

recipe_obj <- recipe(intenzita ~ ., training(d_c_data_splits)) %>%
  step_timeseries_signature(cas_mereni) %>%
  step_rm(cas_mereni) %>%
  step_zv(all_predictors()) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE)

recipe_obj_with_date <- recipe(intenzita ~ ., training(d_c_data_splits)) %>%
  step_timeseries_signature(cas_mereni) %>%
  step_zv(all_predictors()) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE)

model_fit_prophet  <- prophet_reg() %>%
   set_engine('prophet') %>% 
  fit(intenzita ~ ., data = training(d_c_data_splits))

model_fit_arima <- arima_reg() %>%
   set_engine(engine = "auto_arima") %>%
    fit(intenzita ~ ., data = training(d_c_data_splits))
  
model_fit_naive <- naive_reg() %>%
  set_engine("snaive") %>%
  fit(intenzita ~ ., data = training(d_c_data_splits))

workflow_prophet_boost <- workflow() %>%
  add_model(
     prophet_boost() %>% set_engine('prophet_xgboost') 
  ) %>%
  add_recipe(recipe_obj_with_date) %>%
  fit(training(d_c_data_splits))

workflow_arima_boost <- workflow() %>%
  add_model(
     arima_boost(non_seasonal_ar = ) %>% set_engine('arima_xgboost') 
  ) %>%
  add_recipe(recipe_obj_with_date) %>%
  fit(training(d_c_data_splits))


model_fit_nnetar <- nnetar_reg() %>%
  set_engine("nnetar") %>%
  fit(intenzita ~ ., data = training(d_c_data_splits))

workflow_xgboost <- workflow() %>%
  add_model(
    boost_tree("regression") %>% set_engine("xgboost")
  ) %>%
  add_recipe(recipe_obj) %>%
  fit(training(d_c_data_splits))

workflow_randforest <- workflow() %>%
  add_model(
    rand_forest("regression") %>% set_engine("randomForest")
  ) %>%
  add_recipe(recipe_obj) %>%
  fit(training(d_c_data_splits))


models_tbl <- modeltime_table(model_fit_arima, model_fit_prophet, model_fit_naive, workflow_prophet_boost, model_fit_nnetar, workflow_xgboost, workflow_randforest)

calibration_tbl <- models_tbl %>% 
  modeltime_calibrate(new_data = testing(d_c_data_splits))

(forecast_tbl <- calibration_tbl %>%
    modeltime_forecast(
        new_data    = testing(d_c_data_splits),
        actual_data = d_c_data,
        keep_data = TRUE
    ))

forecast_tbl %>% plot_modeltime_forecast()

calibration_tbl %>%
    modeltime_accuracy() %>%
    table_modeltime_accuracy(
        .interactive = FALSE
    )
```


```{r}
models_tbl <- modeltime_table(model_fit_arima, model_fit_prophet, model_fit_naive, workflow_prophet_boost, model_fit_nnetar, workflow_xgboost, workflow_randforest)

calibration_tbl <- models_tbl %>% 
  modeltime_calibrate(new_data = testing(d_c_data_splits$splits[[1]]))

(forecast_tbl <- calibration_tbl %>%
    modeltime_forecast(
        new_data    = testing(d_c_data_splits$splits[[1]]),
        actual_data = d_c_data,
        keep_data = TRUE #Includes the new data (and actual data) as extra columns with the results of the model forecasts
    ))

forecast_tbl %>% plot_modeltime_forecast()

resample_results <- models_tbl %>%
  modeltime_fit_resamples(
    resamples = d_c_data_splits,
    control   = control_resamples(verbose = FALSE)
  )

resample_results %>%
  modeltime_resample_accuracy(summary_fns = mean) %>%
    table_modeltime_accuracy(.interactive = FALSE)
```

```{r}
forecast_tbl %>% 
  filter(str_detect(.model_desc, "ARIMA") | .model_desc == "ACTUAL") %>%
  ggplot(aes(x = cas_mereni, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2020-12-04 03:00","2020-12-04 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "ARIMA")) +
    theme(legend.position="bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model ARIMA", subtitle = "Pro data ze stacionÃ¡rnÃ­ch detektorÅ¯")

ggsave("str_dobr_arima_model.svg")
```

```{r}
forecast_tbl %>% 
  filter(.model_desc == "PROPHET" | .model_desc == "ACTUAL") %>%
  ggplot(aes(x = cas_mereni, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2020-12-04 03:00","2020-12-04 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "Prophet")) +
    theme(legend.position="bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model Prophet", subtitle = "Pro data ze stacionÃ¡rnÃ­ch detektorÅ¯")

ggsave("str_dobr_prophet_model.svg")
```

```{r}
forecast_tbl %>% 
  filter(str_detect(.model_desc, "(SNAIVE)|(ACTUAL)")) %>%
  ggplot(aes(x = cas_mereni, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2020-12-04 03:00","2020-12-04 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "sNaive")) +
    theme(legend.position="bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model sNaive", subtitle = "Pro data ze stacionÃ¡rnÃ­ch detektorÅ¯")

ggsave("str_dobr_snaive_model.svg")
```

```{r}
forecast_tbl %>% 
  filter(str_detect(.model_desc, "(PROPHET W/ XGBOOST ERRORS)|(ACTUAL)")) %>%
  ggplot(aes(x = cas_mereni, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2020-12-04 03:00","2020-12-04 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "Prophet s XGBoost")) +
    theme(legend.position="bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model Prophet s XGBoost", subtitle = "Pro data ze stacionÃ¡rnÃ­ch detektorÅ¯")

ggsave("str_dobr_prophet_boosted_model.svg")
```

```{r}
forecast_tbl %>% 
  filter(str_detect(.model_desc, "(NNAR)|(ACTUAL)")) %>%
  ggplot(aes(x = cas_mereni, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2020-12-04 03:00","2020-12-04 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "NNAR")) +
    theme(legend.position="bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model NNAR", subtitle = "Pro data ze stacionÃ¡rnÃ­ch detektorÅ¯")

ggsave("str_dobr_nnar_model.svg")
```

```{r}
forecast_tbl %>% 
  filter(str_detect(.model_desc, "(RANDOMFOREST)|(ACTUAL)")) %>%
  ggplot(aes(x = cas_mereni, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2020-12-04 03:00","2020-12-04 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "Random Forest")) +
    theme(legend.position="bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model Random Forest", subtitle = "Pro data ze stacionÃ¡rnÃ­ch detektorÅ¯")

ggsave("str_dobr_randomforest_model.svg")
```

```{r}
forecast_tbl %>% 
  filter(.model_desc == "XGBOOST" | .model_desc == "ACTUAL") %>%
  ggplot(aes(x = cas_mereni, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2020-12-04 03:00","2020-12-04 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "XGBoost")) +
    theme(legend.position="bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model XGBoost", subtitle = "Pro data ze stacionÃ¡rnÃ­ch detektorÅ¯")

ggsave("str_dobr_xgboost_model.svg")
```

```{r}
calibration_tbl %>%
    modeltime_accuracy() %>%
    table_modeltime_accuracy(
        .interactive = FALSE
    )
```



### CV

```{r}
d_c_data_splits <- d_c_data %>%
  select(-prumerna_rychlost) %>%
  time_series_cv(date_var = cas_mereni, assess = "1 day", initial = "41 days", skip = "1 days", slice_limit = 7)

recipe_obj <- recipe(intenzita ~ ., training(d_c_data_splits$splits[[1]])) %>%
  step_timeseries_signature(cas_mereni) %>%
  step_rm(cas_mereni) %>%
  step_zv(all_predictors()) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE)

recipe_obj_with_date <- recipe(intenzita ~ ., training(d_c_data_splits$splits[[1]])) %>%
  step_timeseries_signature(cas_mereni) %>%
  step_zv(all_predictors()) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE)

model_fit_prophet  <- prophet_reg() %>%
   set_engine('prophet') %>% 
  fit(intenzita ~ ., data = training(d_c_data_splits$splits[[1]]))


  
model_fit_naive <- naive_reg() %>%
  set_engine("snaive") %>%
  fit(intenzita ~ ., data = training(d_c_data_splits$splits[[1]]))

workflow_prophet_boost <- workflow() %>%
  add_model(
     prophet_boost() %>% set_engine('prophet_xgboost') 
  ) %>%
  add_recipe(recipe_obj_with_date) %>%
  fit(training(d_c_data_splits$splits[[1]]))



model_fit_nnetar <- nnetar_reg() %>%
  set_engine("nnetar") %>%
  fit(intenzita ~ ., data = training(d_c_data_splits$splits[[1]]))

workflow_xgboost <- workflow() %>%
  add_model(
    boost_tree("regression") %>% set_engine("xgboost")
  ) %>%
  add_recipe(recipe_obj) %>%
  fit(training(d_c_data_splits$splits[[1]]))

workflow_randforest <- workflow() %>%
  add_model(
    rand_forest("regression") %>% set_engine("randomForest")
  ) %>%
  add_recipe(recipe_obj) %>%
  fit(training(d_c_data_splits$splits[[1]]))


```

```{r}
d_c_data_splits <- d_c_data %>%
  select(-prumerna_rychlost) %>%
  time_series_cv(date_var = cas_mereni, assess = "1 day", initial = "51 days", skip = "1 days", slice_limit = 5)
```

```{r}
d_c_data_splits <- d_c_data %>%
  select(-prumerna_rychlost) %>%
  filter(cas_mereni < '2020-11-30 00:00:00 ') %>%
  time_series_split(date_var = cas_mereni, assess = "1 day", initial = "60 days", cumulative = FALSE)

d_c_data_splits %>%
tk_time_series_cv_plan() %>%
  plot_time_series_cv_plan(cas_mereni, intenzita, .interactive = TRUE)
```

```{r}
model_fit_arima <- arima_reg() %>%
   set_engine(engine = "auto_arima") %>%
    fit(intenzita ~ ., data = training(d_c_data_splits))

arima_tbl <- modeltime_table(model_fit_arima)

arima_calibration_tbl <- arima_tbl %>% 
  modeltime_calibrate(new_data = testing(d_c_data_splits))

(forecast_tbl <- arima_calibration_tbl %>%
    modeltime_forecast(
        new_data    = testing(d_c_data_splits),
        actual_data = d_c_data,
        keep_data = TRUE 
    ))

forecast_tbl %>% plot_modeltime_forecast()

arima_calibration_tbl %>%
    modeltime_accuracy() %>%
    table_modeltime_accuracy(
        .interactive = FALSE
    )
```

```{r}
arima_tbl <- modeltime_table(model_fit_arima)

arima_calibration_tbl <- arima_tbl %>% 
  modeltime_calibrate(new_data = testing(d_c_data_splits$splits[[1]]))

(forecast_tbl <- arima_calibration_tbl %>%
    modeltime_forecast(
        new_data    = testing(d_c_data_splits$splits[[1]]),
        actual_data = d_c_data,
        keep_data = TRUE 
    ))

forecast_tbl %>% plot_modeltime_forecast()

resample_results <- arima_tbl %>%
  modeltime_fit_resamples(
    resamples = d_c_data_splits,
    control   = control_resamples(verbose = FALSE)
  )

resample_results %>%
  modeltime_resample_accuracy(summary_fns = mean) %>%
    table_modeltime_accuracy(.interactive = FALSE)
```


```{r}
models_tbl <- modeltime_table(model_fit_prophet, model_fit_naive, workflow_prophet_boost, model_fit_nnetar, workflow_xgboost, workflow_randforest)

calibration_tbl <- models_tbl %>% 
  modeltime_calibrate(new_data = testing(d_c_data_splits$splits[[1]]))

(forecast_tbl <- calibration_tbl %>%
    modeltime_forecast(
        new_data    = testing(d_c_data_splits$splits[[1]]),
        actual_data = d_c_data,
        keep_data = TRUE
    ))

forecast_tbl %>% plot_modeltime_forecast()

resample_results <- models_tbl %>%
  modeltime_fit_resamples(
    resamples = d_c_data_splits,
    control   = control_resamples(verbose = FALSE)
  )

resample_results %>%
  modeltime_resample_accuracy(summary_fns = mean) %>%
    table_modeltime_accuracy(.interactive = FALSE)
```

### KNN

```{r}
d_c_data_knn <- d_c_data %>%
  filter(cas_mereni > "2020-09-03 12:55:00" & cas_mereni <= "2020-09-03 22:55:00")

d_c_knn_split <- "2020-09-03 17:55:00"

d_c_data_knn_training <- d_c_data_knn %>%
  filter(cas_mereni <= d_c_knn_split ) %>%
  select(c(cas_mereni, intenzita)) %>%
  read.zoo() %>%
  as.ts()

d_c_data_knn_testing <- d_c_data_knn %>%
  filter(cas_mereni > d_c_knn_split) %>%
  select(c(cas_mereni, intenzita)) %>%
  mutate(type = "testing")

knn_model <- knn_forecasting(timeS = d_c_data_knn_training, h = 6, lags = 1:6, k = 7)

knn_plot <- autoplot(knn_model)

knn_pred_tbl <- knn_model$prediction %>%
  tk_tbl(rename_index = "cas_mereni") %>%
  mutate(cas_mereni = as_datetime(cas_mereni)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

d_c_data_knn %>%
  select(c(cas_mereni, intenzita)) %>%
  mutate(type = "data") %>%
  bind_rows(knn_pred_tbl) %>%
ggplot(aes(x = cas_mereni, y = intenzita, color = type)) +
  geom_line() +
  theme_classic() +
  scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "K-NN")) +
  scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2020-09-03 13:00:002","2020-09-03 18:30:00"), format = "%Y-%m-%d %H:%M"))) +
  labs(x = "Äas", y = "DopravnÃ­ intenzita [voz/hod]", color = "", title = "Model K-NN")
  

```
### KNN manualni cv

```{r}
d_c_knn_begin <- "2020-09-07 00:00:00"
d_c_knn_end <- "2020-11-30 23:59:00"
d_c_knn_split <- "2020-11-30 00:00:00"

d_c_data_knn <- d_c_data %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni <= d_c_knn_end)



d_c_data_knn_training <- d_c_data_knn %>%
  filter(cas_mereni <= d_c_knn_split ) %>%
  select(c(cas_mereni, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

d_c_data_knn_testing <- d_c_data %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(cas_mereni, intenzita)) %>%
  mutate(type = "testing")

knn_model <- knn_forecasting(timeS = d_c_data_knn_training, h = 260, lags = 1:260, k = 7)

knn_pred_tbl <- knn_model$prediction %>%
  tk_tbl(rename_index = "cas_mereni") %>%
  mutate(cas_mereni = as_datetime(cas_mereni)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

d_c_knn_eval <- d_c_data_knn_testing %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni < d_c_knn_end) %>%
  select(-type) %>%
  left_join(knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "cas_mereni") 

knn_mae_dobr <- NA
knn_mape_dobr <- NA
knn_mase_dobr <- NA
knn_smape_dobr <- NA
knn_rmse_dobr <- NA
knn_r2_dobr <- NA

knn_mae_dobr %<>% append((d_c_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_dobr %<>% append((d_c_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_dobr %<>% append((d_c_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_dobr %<>% append((d_c_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_dobr %<>% append((d_c_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_dobr %<>% append((d_c_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

#--------
d_c_knn_begin <- "2020-09-06 00:00:00"
d_c_knn_end <- "2020-11-29 23:59:00"
d_c_knn_split <- "2020-11-29 00:00:00"

d_c_data_knn <- d_c_data %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni <= d_c_knn_end)



d_c_data_knn_training <- d_c_data_knn %>%
  filter(cas_mereni <= d_c_knn_split ) %>%
  select(c(cas_mereni, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

d_c_data_knn_testing <- d_c_data %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(cas_mereni, intenzita)) %>%
  mutate(type = "testing")

knn_model <- knn_forecasting(timeS = d_c_data_knn_training, h = 260, lags = 1:260, k = 7)

knn_pred_tbl <- knn_model$prediction %>%
  tk_tbl(rename_index = "cas_mereni") %>%
  mutate(cas_mereni = as_datetime(cas_mereni)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

d_c_knn_eval <- d_c_data_knn_testing %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni < d_c_knn_end) %>%
  select(-type) %>%
  left_join(knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "cas_mereni") 


knn_mae_dobr %<>% append((d_c_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_dobr %<>% append((d_c_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_dobr %<>% append((d_c_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_dobr %<>% append((d_c_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_dobr %<>% append((d_c_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_dobr %<>% append((d_c_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

#--------
d_c_knn_begin <- "2020-09-05 00:00:00"
d_c_knn_end <- "2020-11-28 23:59:00"
d_c_knn_split <- "2020-11-28 00:00:00"

d_c_data_knn <- d_c_data %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni <= d_c_knn_end)



d_c_data_knn_training <- d_c_data_knn %>%
  filter(cas_mereni <= d_c_knn_split ) %>%
  select(c(cas_mereni, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

d_c_data_knn_testing <- d_c_data %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(cas_mereni, intenzita)) %>%
  mutate(type = "testing")

knn_model <- knn_forecasting(timeS = d_c_data_knn_training, h = 260, lags = 1:260, k = 7)

knn_pred_tbl <- knn_model$prediction %>%
  tk_tbl(rename_index = "cas_mereni") %>%
  mutate(cas_mereni = as_datetime(cas_mereni)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

d_c_knn_eval <- d_c_data_knn_testing %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni < d_c_knn_end) %>%
  select(-type) %>%
  left_join(knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "cas_mereni") 


knn_mae_dobr %<>% append((d_c_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_dobr %<>% append((d_c_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_dobr %<>% append((d_c_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_dobr %<>% append((d_c_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_dobr %<>% append((d_c_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_dobr %<>% append((d_c_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

#--------
d_c_knn_begin <- "2020-09-04 00:00:00"
d_c_knn_end <- "2020-11-27 23:59:00"
d_c_knn_split <- "2020-11-27 00:00:00"

d_c_data_knn <- d_c_data %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni <= d_c_knn_end)



d_c_data_knn_training <- d_c_data_knn %>%
  filter(cas_mereni <= d_c_knn_split ) %>%
  select(c(cas_mereni, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

d_c_data_knn_testing <- d_c_data %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(cas_mereni, intenzita)) %>%
  mutate(type = "testing")

knn_model <- knn_forecasting(timeS = d_c_data_knn_training, h = 260, lags = 1:260, k = 7)

knn_pred_tbl <- knn_model$prediction %>%
  tk_tbl(rename_index = "cas_mereni") %>%
  mutate(cas_mereni = as_datetime(cas_mereni)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

d_c_knn_eval <- d_c_data_knn_testing %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni < d_c_knn_end) %>%
  select(-type) %>%
  left_join(knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "cas_mereni") 


knn_mae_dobr %<>% append((d_c_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_dobr %<>% append((d_c_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_dobr %<>% append((d_c_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_dobr %<>% append((d_c_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_dobr %<>% append((d_c_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_dobr %<>% append((d_c_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

#--------
d_c_knn_begin <- "2020-09-03 00:00:00"
d_c_knn_end <- "2020-11-26 23:59:00"
d_c_knn_split <- "2020-11-26 00:00:00"

d_c_data_knn <- d_c_data %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni <= d_c_knn_end)



d_c_data_knn_training <- d_c_data_knn %>%
  filter(cas_mereni <= d_c_knn_split ) %>%
  select(c(cas_mereni, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

d_c_data_knn_testing <- d_c_data %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(cas_mereni, intenzita)) %>%
  mutate(type = "testing")

knn_model <- knn_forecasting(timeS = d_c_data_knn_training, h = 260, lags = 1:260, k = 7)

knn_pred_tbl <- knn_model$prediction %>%
  tk_tbl(rename_index = "cas_mereni") %>%
  mutate(cas_mereni = as_datetime(cas_mereni)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

d_c_knn_eval <- d_c_data_knn_testing %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni < d_c_knn_end) %>%
  select(-type) %>%
  left_join(knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "cas_mereni") 


knn_mae_dobr %<>% append((d_c_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_dobr %<>% append((d_c_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_dobr %<>% append((d_c_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_dobr %<>% append((d_c_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_dobr %<>% append((d_c_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_dobr %<>% append((d_c_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

#--------
d_c_knn_begin <- "2020-09-02 00:00:00"
d_c_knn_end <- "2020-11-25 23:59:00"
d_c_knn_split <- "2020-11-25 00:00:00"

d_c_data_knn <- d_c_data %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni <= d_c_knn_end)



d_c_data_knn_training <- d_c_data_knn %>%
  filter(cas_mereni <= d_c_knn_split ) %>%
  select(c(cas_mereni, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

d_c_data_knn_testing <- d_c_data %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(cas_mereni, intenzita)) %>%
  mutate(type = "testing")

knn_model <- knn_forecasting(timeS = d_c_data_knn_training, h = 260, lags = 1:260, k = 7)

knn_pred_tbl <- knn_model$prediction %>%
  tk_tbl(rename_index = "cas_mereni") %>%
  mutate(cas_mereni = as_datetime(cas_mereni)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

d_c_knn_eval <- d_c_data_knn_testing %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni < d_c_knn_end) %>%
  select(-type) %>%
  left_join(knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "cas_mereni") 


knn_mae_dobr %<>% append((d_c_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_dobr %<>% append((d_c_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_dobr %<>% append((d_c_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_dobr %<>% append((d_c_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_dobr %<>% append((d_c_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_dobr %<>% append((d_c_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

#--------
d_c_knn_begin <- "2020-09-01 00:00:00"
d_c_knn_end <- "2020-11-24 23:59:00"
d_c_knn_split <- "2020-11-24 00:00:00"

d_c_data_knn <- d_c_data %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni <= d_c_knn_end)



d_c_data_knn_training <- d_c_data_knn %>%
  filter(cas_mereni <= d_c_knn_split ) %>%
  select(c(cas_mereni, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

d_c_data_knn_testing <- d_c_data %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(cas_mereni, intenzita)) %>%
  mutate(type = "testing")

knn_model <- knn_forecasting(timeS = d_c_data_knn_training, h = 260, lags = 1:260, k = 7)

knn_pred_tbl <- knn_model$prediction %>%
  tk_tbl(rename_index = "cas_mereni") %>%
  mutate(cas_mereni = as_datetime(cas_mereni)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

d_c_knn_eval <- d_c_data_knn_testing %>%
  filter(cas_mereni > d_c_knn_begin & cas_mereni < d_c_knn_end) %>%
  select(-type) %>%
  left_join(knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "cas_mereni") 


knn_mae_dobr %<>% append((d_c_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_dobr %<>% append((d_c_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_dobr %<>% append((d_c_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_dobr %<>% append((d_c_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_dobr %<>% append((d_c_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_dobr %<>% append((d_c_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

```

```{r}
knn_mae_dobr %<>% na.omit() %>% mean()
knn_mape_dobr %<>% na.omit() %>% mean()
knn_mase_dobr %<>% na.omit() %>% mean()
knn_smape_dobr %<>% na.omit() %>% mean()
knn_rmse_dobr %<>% na.omit() %>% mean()
knn_r2_dobr %<>% na.omit() %>% mean()

cat("MAE = ", knn_mae_dobr, "\n")
cat("MAPE = ", knn_mape_dobr, "\n")
cat("MASE = ", knn_mase_dobr, "\n")
cat("SMAPE = ", knn_smape_dobr, "\n")
cat("RMSE = ", knn_rmse_dobr, "\n")
cat("R^2 = ", knn_r2_dobr)

## FCD Predikce 

```{r}
#priprava pro cross validaci
fcd_data_rounded_splits <- fcd_data_rounded %>%
  time_series_split(date_var = datum, assess = "1 day", cumulative = TRUE)
  #time_series_cv(date_var = cas_mereni, assess = "1 day", initial = "8 days", skip = "1 days", slice_limit = 7) 
  
```

### FCD predikce bez cross validace

```{r}
fcd_data_rounded_splits <- fcd_data_rounded %>%
  #filter(datum < '2019-11-29 00:00:00') %>%
  time_series_split(date_var = datum, assess = "1 days", cumulative = TRUE)

recipe_obj <- recipe(intenzita ~ ., training(fcd_data_rounded_splits)) %>%
  step_timeseries_signature(datum) %>%
  step_rm(datum) %>%
  step_zv(all_predictors()) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE)

recipe_obj_with_date <- recipe(intenzita ~ ., training(fcd_data_rounded_splits)) %>%
  step_timeseries_signature(datum) %>%
  step_zv(all_predictors()) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE)


model_fit_prophet  <- prophet_reg() %>%
    set_engine('prophet') %>% 
  fit(intenzita ~ ., data = training(fcd_data_rounded_splits))

workflow_prophet_boost <- workflow() %>%
  add_model(
     prophet_boost() %>% set_engine('prophet_xgboost') 
  ) %>%
  add_recipe(recipe_obj_with_date) %>%
  fit(training(fcd_data_rounded_splits))

model_fit_arima <- arima_reg() %>%
    set_engine(engine = "auto_arima") %>%
    fit(intenzita ~ ., data = training(fcd_data_rounded_splits))

model_fit_naive <- naive_reg() %>%
  set_engine("snaive") %>%
  fit(intenzita ~ datum, data = training(fcd_data_rounded_splits))

model_fit_nnetar <- nnetar_reg() %>%
  set_engine("nnetar") %>%
  fit(intenzita ~ datum, data = training(fcd_data_rounded_splits))

workflow_xgboost <- workflow() %>%
  add_model(
    boost_tree("regression") %>% set_engine("xgboost")
  ) %>%
  add_recipe(recipe_obj) %>%
  fit(training(fcd_data_rounded_splits))

workflow_randforest <- workflow() %>%
  add_model(
    rand_forest("regression") %>% set_engine("randomForest")
  ) %>%
  add_recipe(recipe_obj) %>%
  fit(training(fcd_data_rounded_splits))

#model_fit_arima_boosted <- arima_boost(
#    min_n = 2,
#    learn_rate = 0.015
#) %>%
#    set_engine(engine = "auto_arima_xgboost") %>%
#    fit(n_vozidel ~ cas_mereni + as.numeric(cas_mereni) + factor(day(cas_mereni)),
#        data = training(d_c_senzor_5_den_splits))

#models_tbl <- modeltime_table(model_fit_prophet, model_fit_arima_no_boost, model_fit_naive, model_fit_nnetar)
models_tbl <- modeltime_table(model_fit_arima, model_fit_prophet, model_fit_naive, model_fit_nnetar, workflow_prophet_boost, workflow_xgboost, workflow_randforest)


calibration_tbl <- models_tbl %>% 
  modeltime_calibrate(new_data = testing(fcd_data_rounded_splits))

(forecast_tbl <- calibration_tbl %>%
    modeltime_forecast(
        new_data    = testing(fcd_data_rounded_splits),
        actual_data = fcd_data_rounded,
        keep_data = TRUE #Includes the new data (and actual data) as extra columns with the results of the model forecasts
    ))

forecast_tbl %>% plot_modeltime_forecast()
```

```{r}
calibration_tbl %>%
    modeltime_accuracy() %>%
    table_modeltime_accuracy(
        .interactive = FALSE
    )
```


```{r}
forecast_tbl %>% 
  filter(str_detect(.model_desc, "(ARIMA)|(ACTUAL)")) %>%
  ggplot(aes(x = datum, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2019-11-30 03:00","2020-11-31 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "ARIMA")) +
    theme(legend.position="bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model ARIMA", subtitle = "Pro FCD z DobÅichovic")

ggsave("fcd_dobr_arima_model.svg")
```

```{r}
forecast_tbl %>% 
  filter(.model_desc == "PROPHET" | .model_desc == "ACTUAL") %>%
  ggplot(aes(x = datum, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2019-11-30 03:00","2020-11-31 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "Prophet")) +
    theme(legend.position="bottom") +
    labs(x = "Äas [hh:mm] ", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model Prophet", subtitle = "Pro FCD z DobÅichovic")

ggsave("fcd_dobr_prophet_model.svg")
```
```{r}
forecast_tbl %>% 
  filter(str_detect(.model_desc, "(SNAIVE)|(ACTUAL)")) %>%
  ggplot(aes(x = datum, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2019-11-30 03:00","2020-11-31 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "sNaive")) +
    theme(legend.position="bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model sNaive", subtitle = "Pro FCD z DobÅichovic")

ggsave("fcd_dobr_snaive_model.svg")
```

```{r}
forecast_tbl %>% 
  filter(.model_desc == "RANDOMFOREST" | .model_desc == "ACTUAL") %>%
  ggplot(aes(x = datum, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2019-11-30 03:00","2020-11-31 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "Random Forest")) +
    theme(legend.position = "bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model Random Forest", subtitle = "Pro FCD z DobÅichovic")

ggsave("fcd_dobr_randomforest_model.svg")
```

```{r}
forecast_tbl %>% 
  filter(str_detect(.model_desc, "(PROPHET W/ XGBOOST ERRORS)|(ACTUAL)")) %>%
  ggplot(aes(x = datum, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2019-11-30 03:00","2020-11-31 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "Prophet s XGBoost")) +
    theme(legend.position="bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model Prophet s XGBoost", subtitle = "Pro FCD z DobÅichovic")

ggsave("fcd_dobr_prophet_boosted_model.svg")
```

```{r}
forecast_tbl %>% 
  filter(.model_desc == "XGBOOST" | .model_desc == "ACTUAL") %>%
  ggplot(aes(x = datum, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2019-11-30 03:00","2020-11-31 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "XGBoost")) +
    theme(legend.position = "bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model XGBoost", subtitle = "Pro FCD z DobÅichovic")

ggsave("fcd_dobr_xgboost_model.svg")
```

```{r}
forecast_tbl %>% 
  filter(str_detect(.model_desc, "(NNAR)|(ACTUAL)")) %>%
  ggplot(aes(x = datum, y = .value)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, color = .model_desc), fill= "gray70",  linetype = "solid", linewidth = 0, alpha = .5, show.legend = FALSE) +
    geom_line(aes(color = .model_desc)) +
    theme_classic() +
    scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2019-11-30 03:00","2020-11-31 23:00"), format = "%Y-%m-%d %H:%M"))) +
    scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "NNAR")) +
    theme(legend.position="bottom") +
    labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model NNAR", subtitle = "Pro FCD z DobÅichovic")

ggsave("fcd_dobr_nnar_model.svg")
```

### FCD Predikce s cross validaci

```{r}
#priprava pro cross validaci
fcd_data_rounded_splits <- fcd_data_rounded %>%
  time_series_cv(date_var = datum, assess = "1 day", initial = "80 days", skip = "1 days", slice_limit = 7) 
  
```

```{r}
recipe_obj <- recipe(intenzita ~ ., training(fcd_data_rounded_splits$splits[[1]])) %>%
  step_timeseries_signature(datum) %>%
  step_rm(datum) %>%
  step_zv(all_predictors()) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE)

recipe_obj_with_date <- recipe(intenzita ~ ., training(fcd_data_rounded_splits$splits[[1]])) %>%
  step_timeseries_signature(datum) %>%
  step_zv(all_predictors()) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE)


model_fit_prophet  <- prophet_reg() %>%
    set_engine('prophet') %>% 
  fit(intenzita ~ ., data = training(fcd_data_rounded_splits$splits[[1]]))

workflow_prophet_boost <- workflow() %>%
  add_model(
     prophet_boost() %>% set_engine('prophet_xgboost') 
  ) %>%
  add_recipe(recipe_obj_with_date) %>%
  fit(training(fcd_data_rounded_splits$splits[[1]]))

model_fit_arima <- arima_reg() %>%
    set_engine(engine = "auto_arima") %>%
    fit(intenzita ~ ., data = training(fcd_data_rounded_splits$splits[[1]]))

model_fit_naive <- naive_reg() %>%
  set_engine("snaive") %>%
  fit(intenzita ~ datum, data = training(fcd_data_rounded_splits$splits[[1]]))

model_fit_nnetar <- nnetar_reg() %>%
  set_engine("nnetar") %>%
  fit(intenzita ~ datum, data = training(fcd_data_rounded_splits$splits[[1]]))

workflow_xgboost <- workflow() %>%
  add_model(
    boost_tree("regression") %>% set_engine("xgboost")
  ) %>%
  add_recipe(recipe_obj) %>%
  fit(training(fcd_data_rounded_splits$splits[[1]]))

workflow_randforest <- workflow() %>%
  add_model(
    rand_forest("regression") %>% set_engine("randomForest")
  ) %>%
  add_recipe(recipe_obj) %>%
  fit(training(fcd_data_rounded_splits$splits[[1]]))

#model_fit_arima_boosted <- arima_boost(
#    min_n = 2,
#    learn_rate = 0.015
#) %>%
#    set_engine(engine = "auto_arima_xgboost") %>%
#    fit(n_vozidel ~ cas_mereni + as.numeric(cas_mereni) + factor(day(cas_mereni)),
#        data = training(d_c_senzor_5_den_splits))

#models_tbl <- modeltime_table(model_fit_prophet, model_fit_arima_no_boost, model_fit_naive, model_fit_nnetar)
```

```{r}
models_tbl <- modeltime_table(model_fit_arima, model_fit_prophet, model_fit_naive, workflow_prophet_boost, model_fit_nnetar, workflow_xgboost, workflow_randforest)

calibration_tbl <- models_tbl %>% 
  modeltime_calibrate(new_data = testing(fcd_data_rounded_splits$splits[[1]]))

(forecast_tbl <- calibration_tbl %>%
    modeltime_forecast(
        new_data    = testing(fcd_data_rounded_splits$splits[[1]]),
        actual_data = fcd_data_rounded,
        keep_data = TRUE #Includes the new data (and actual data) as extra columns with the results of the model forecasts
    ))

resample_results <- models_tbl %>%
  modeltime_fit_resamples(
    resamples = fcd_data_rounded_splits,
    control   = control_resamples(verbose = FALSE)
  )
```

```{r}
forecast_tbl %>% plot_modeltime_forecast()

resample_results %>%
  modeltime_resample_accuracy(summary_fns = mean) %>%
    table_modeltime_accuracy(.interactive = FALSE)
```



### FCD KNN predikce

```{r}
fcd_data_knn <- fcd_data_rounded %>%
  filter(datum > "2019-09-01 00:00:00" & datum <= "2019-11-30 23:59:00")

fcd_data_knn_split <- "2019-11-30 11:00:00"

fcd_data_knn_training <- fcd_data_knn %>%
  filter(datum <= fcd_data_knn_split ) %>%
  select(c(datum, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

fcd_data_knn_testing <- fcd_data_rounded %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(datum, intenzita)) %>%
  mutate(type = "testing")

fcd_knn_model <- knn_forecasting(timeS = fcd_data_knn_training, h = 260, lags = 1:260, k = 7)

fcd_knn_plot <- autoplot(fcd_knn_model)

fcd_knn_pred_tbl <- fcd_knn_model$prediction %>%
  tk_tbl(rename_index = "datum") %>%
  mutate(datum = as_datetime(datum)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

fcd_data_knn_testing %>%
  mutate(type = "data") %>%
  bind_rows(fcd_knn_pred_tbl) %>%
ggplot(aes(x = datum, y = intenzita, color = type)) +
  geom_line() +
  theme_classic() +
  scale_color_aaas(labels = c("SkuteÄnÃ¡ data", "K-NN")) +
  scale_x_datetime(date_labels = "%H:%M", limits = as.POSIXct(strptime(c("2019-11-30 00:00:00","2019-11-30 23:59:00"), format = "%Y-%m-%d %H:%M"))) +
  labs(x = "Äas [hh:mm]", y = "DopravnÃ­ intenzita [voz/h]", color = "", title = "Model K-NN", subtitle = "Pro FCD z DobÅichovic") +
  theme(legend.position = "bottom")

ggsave("fcd_dobr_knn.svg")
  
fcd_knn_plot

#d_c_data %>%
#  filter(cas_mereni > "2020-09-03 12:55:00" & cas_mereni <= "2020-09-03 22:55:00") %>%
#  select(c(cas_mereni, intenzita)) %>%
#  mutate(type = "training") %>%
#  bind_rows(knn_pred_tbl) %>%
#ggplot(aes(x = cas_mereni, y = intenzita, color = type)) +
#  geom_line()

```

```{r}
fcd_knn_eval <- fcd_data_knn_testing %>%
  filter(datum > '2019-11-30 00:00:00' & datum < '2019-11-30 23:59:00') %>%
  select(-type) %>%
  left_join(fcd_knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "datum") 

fcd_knn_eval %>% yardstick::mae(intenzita, preds)

fcd_knn_eval %>% yardstick::mape(intenzita, preds)

fcd_knn_eval %>% yardstick::mase(intenzita, preds)

fcd_knn_eval %>% yardstick::smape(intenzita, preds)

fcd_knn_eval %>% yardstick::rmse(intenzita, preds)

fcd_knn_eval %>% yardstick::rsq(intenzita, preds)
```

### FCD KNN manualni cv

```{r}
fcd_knn_begin <- "2019-09-07 00:00:00"
fcd_knn_end <- "2019-11-30 23:59:00"
fcd_data_knn_split <- "2019-11-30 00:00:00"

fcd_data_knn <- fcd_data_rounded %>%
  filter(datum > fcd_knn_begin & datum <= fcd_knn_end)


fcd_data_knn_training <- fcd_data_knn %>%
  filter(datum <= fcd_data_knn_split ) %>%
  select(c(datum, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

fcd_data_knn_testing <- fcd_data_rounded %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(datum, intenzita)) %>%
  mutate(type = "testing")

fcd_knn_model <- knn_forecasting(timeS = fcd_data_knn_training, h = 260, lags = 1:260, k = 7)

fcd_knn_plot <- autoplot(fcd_knn_model)

fcd_knn_pred_tbl <- fcd_knn_model$prediction %>%
  tk_tbl(rename_index = "datum") %>%
  mutate(datum = as_datetime(datum)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

fcd_knn_eval <- fcd_data_knn_testing %>%
  filter(datum > fcd_knn_begin & datum < fcd_knn_end) %>%
  select(-type) %>%
  left_join(fcd_knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "datum") 

knn_mae_fcd <- NA
knn_mape_fcd <- NA
knn_mase_fcd <- NA
knn_smape_fcd <- NA
knn_rmse_fcd <- NA
knn_r2_fcd <- NA

knn_mae_fcd %<>% append((fcd_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_fcd %<>% append((fcd_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_fcd %<>% append((fcd_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_fcd %<>% append((fcd_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_fcd %<>% append((fcd_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_fcd %<>% append((fcd_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

#----

fcd_knn_begin <- "2019-09-06 00:00:00"
fcd_knn_end <- "2019-11-29 23:59:00"
fcd_data_knn_split <- "2019-11-29 00:00:00"

fcd_data_knn <- fcd_data_rounded %>%
  filter(datum > fcd_knn_begin & datum <= fcd_knn_end)


fcd_data_knn_training <- fcd_data_knn %>%
  filter(datum <= fcd_data_knn_split ) %>%
  select(c(datum, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

fcd_data_knn_testing <- fcd_data_rounded %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(datum, intenzita)) %>%
  mutate(type = "testing")

fcd_knn_model <- knn_forecasting(timeS = fcd_data_knn_training, h = 260, lags = 1:260, k = 7)

fcd_knn_plot <- autoplot(fcd_knn_model)

fcd_knn_pred_tbl <- fcd_knn_model$prediction %>%
  tk_tbl(rename_index = "datum") %>%
  mutate(datum = as_datetime(datum)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

fcd_knn_eval <- fcd_data_knn_testing %>%
  filter(datum > fcd_knn_begin & datum < fcd_knn_end) %>%
  select(-type) %>%
  left_join(fcd_knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "datum") 

knn_mae_fcd %<>% append((fcd_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_fcd %<>% append((fcd_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_fcd %<>% append((fcd_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_fcd %<>% append((fcd_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_fcd %<>% append((fcd_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_fcd %<>% append((fcd_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

#----

fcd_knn_begin <- "2019-09-05 00:00:00"
fcd_knn_end <- "2019-11-28 23:59:00"
fcd_data_knn_split <- "2019-11-28 00:00:00"

fcd_data_knn <- fcd_data_rounded %>%
  filter(datum > fcd_knn_begin & datum <= fcd_knn_end)


fcd_data_knn_training <- fcd_data_knn %>%
  filter(datum <= fcd_data_knn_split ) %>%
  select(c(datum, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

fcd_data_knn_testing <- fcd_data_rounded %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(datum, intenzita)) %>%
  mutate(type = "testing")

fcd_knn_model <- knn_forecasting(timeS = fcd_data_knn_training, h = 260, lags = 1:260, k = 7)

fcd_knn_plot <- autoplot(fcd_knn_model)

fcd_knn_pred_tbl <- fcd_knn_model$prediction %>%
  tk_tbl(rename_index = "datum") %>%
  mutate(datum = as_datetime(datum)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

fcd_knn_eval <- fcd_data_knn_testing %>%
  filter(datum > fcd_knn_begin & datum < fcd_knn_end) %>%
  select(-type) %>%
  left_join(fcd_knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "datum") 

knn_mae_fcd %<>% append((fcd_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_fcd %<>% append((fcd_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_fcd %<>% append((fcd_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_fcd %<>% append((fcd_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_fcd %<>% append((fcd_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_fcd %<>% append((fcd_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

#----

fcd_knn_begin <- "2019-09-04 00:00:00"
fcd_knn_end <- "2019-11-27 23:59:00"
fcd_data_knn_split <- "2019-11-27 00:00:00"

fcd_data_knn <- fcd_data_rounded %>%
  filter(datum > fcd_knn_begin & datum <= fcd_knn_end)


fcd_data_knn_training <- fcd_data_knn %>%
  filter(datum <= fcd_data_knn_split ) %>%
  select(c(datum, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

fcd_data_knn_testing <- fcd_data_rounded %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(datum, intenzita)) %>%
  mutate(type = "testing")

fcd_knn_model <- knn_forecasting(timeS = fcd_data_knn_training, h = 260, lags = 1:260, k = 7)

fcd_knn_plot <- autoplot(fcd_knn_model)

fcd_knn_pred_tbl <- fcd_knn_model$prediction %>%
  tk_tbl(rename_index = "datum") %>%
  mutate(datum = as_datetime(datum)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

fcd_knn_eval <- fcd_data_knn_testing %>%
  filter(datum > fcd_knn_begin & datum < fcd_knn_end) %>%
  select(-type) %>%
  left_join(fcd_knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "datum") 

knn_mae_fcd %<>% append((fcd_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_fcd %<>% append((fcd_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_fcd %<>% append((fcd_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_fcd %<>% append((fcd_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_fcd %<>% append((fcd_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_fcd %<>% append((fcd_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

#----

fcd_knn_begin <- "2019-09-03 00:00:00"
fcd_knn_end <- "2019-11-27 23:59:00"
fcd_data_knn_split <- "2019-11-27 00:00:00"

fcd_data_knn <- fcd_data_rounded %>%
  filter(datum > fcd_knn_begin & datum <= fcd_knn_end)


fcd_data_knn_training <- fcd_data_knn %>%
  filter(datum <= fcd_data_knn_split ) %>%
  select(c(datum, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

fcd_data_knn_testing <- fcd_data_rounded %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(datum, intenzita)) %>%
  mutate(type = "testing")

fcd_knn_model <- knn_forecasting(timeS = fcd_data_knn_training, h = 260, lags = 1:260, k = 7)

fcd_knn_plot <- autoplot(fcd_knn_model)

fcd_knn_pred_tbl <- fcd_knn_model$prediction %>%
  tk_tbl(rename_index = "datum") %>%
  mutate(datum = as_datetime(datum)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

fcd_knn_eval <- fcd_data_knn_testing %>%
  filter(datum > fcd_knn_begin & datum < fcd_knn_end) %>%
  select(-type) %>%
  left_join(fcd_knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "datum") 

knn_mae_fcd %<>% append((fcd_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_fcd %<>% append((fcd_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_fcd %<>% append((fcd_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_fcd %<>% append((fcd_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_fcd %<>% append((fcd_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_fcd %<>% append((fcd_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

#----

fcd_knn_begin <- "2019-09-02 00:00:00"
fcd_knn_end <- "2019-11-27 23:59:00"
fcd_data_knn_split <- "2019-11-27 00:00:00"

fcd_data_knn <- fcd_data_rounded %>%
  filter(datum > fcd_knn_begin & datum <= fcd_knn_end)


fcd_data_knn_training <- fcd_data_knn %>%
  filter(datum <= fcd_data_knn_split ) %>%
  select(c(datum, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

fcd_data_knn_testing <- fcd_data_rounded %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(datum, intenzita)) %>%
  mutate(type = "testing")

fcd_knn_model <- knn_forecasting(timeS = fcd_data_knn_training, h = 260, lags = 1:260, k = 7)

fcd_knn_plot <- autoplot(fcd_knn_model)

fcd_knn_pred_tbl <- fcd_knn_model$prediction %>%
  tk_tbl(rename_index = "datum") %>%
  mutate(datum = as_datetime(datum)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

fcd_knn_eval <- fcd_data_knn_testing %>%
  filter(datum > fcd_knn_begin & datum < fcd_knn_end) %>%
  select(-type) %>%
  left_join(fcd_knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "datum") 

knn_mae_fcd %<>% append((fcd_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_fcd %<>% append((fcd_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_fcd %<>% append((fcd_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_fcd %<>% append((fcd_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_fcd %<>% append((fcd_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_fcd %<>% append((fcd_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

#----

fcd_knn_begin <- "2019-09-01 00:00:00"
fcd_knn_end <- "2019-11-26 23:59:00"
fcd_data_knn_split <- "2019-11-26 00:00:00"

fcd_data_knn <- fcd_data_rounded %>%
  filter(datum > fcd_knn_begin & datum <= fcd_knn_end)


fcd_data_knn_training <- fcd_data_knn %>%
  filter(datum <= fcd_data_knn_split ) %>%
  select(c(datum, intenzita)) %>%
  read.zoo() %>%
  as.ts() %>%
  zoo() %>%
  na.fill(0) %>%
  as.ts()

fcd_data_knn_testing <- fcd_data_rounded %>%
  #filter(cas_mereni > d_c_knn_split) %>%
  select(c(datum, intenzita)) %>%
  mutate(type = "testing")

fcd_knn_model <- knn_forecasting(timeS = fcd_data_knn_training, h = 260, lags = 1:260, k = 7)

fcd_knn_plot <- autoplot(fcd_knn_model)

fcd_knn_pred_tbl <- fcd_knn_model$prediction %>%
  tk_tbl(rename_index = "datum") %>%
  mutate(datum = as_datetime(datum)) %>%
  rename(intenzita = value) %>%
  mutate(type = "pred")

fcd_knn_eval <- fcd_data_knn_testing %>%
  filter(datum > fcd_knn_begin & datum < fcd_knn_end) %>%
  select(-type) %>%
  left_join(fcd_knn_pred_tbl %>% rename(preds = intenzita) %>% select(-type), by = "datum") 

knn_mae_fcd %<>% append((fcd_knn_eval %>% yardstick::mae(intenzita, preds))$.estimate) 

knn_mape_fcd %<>% append((fcd_knn_eval %>% yardstick::mape(intenzita, preds))$.estimate)

knn_mase_fcd %<>% append((fcd_knn_eval %>% yardstick::mase(intenzita, preds))$.estimate)

knn_smape_fcd %<>% append((fcd_knn_eval %>% yardstick::smape(intenzita, preds))$.estimate)

knn_rmse_fcd %<>% append((fcd_knn_eval %>% yardstick::rmse(intenzita, preds))$.estimate)

knn_r2_fcd %<>% append((fcd_knn_eval %>% yardstick::rsq(intenzita, preds))$.estimate)

```

```{r}
knn_mae_fcd %<>% na.omit() %>% mean()
knn_mape_fcd %<>% na.omit() %>% mean()
knn_mase_fcd %<>% na.omit() %>% mean()
knn_smape_fcd %<>% na.omit() %>% mean()
knn_rmse_fcd %<>% na.omit() %>% mean()
knn_r2_fcd %<>% na.omit() %>% mean()

cat("MAE = ", knn_mae_fcd, "\n")
cat("MAPE = ", knn_mape_fcd, "\n")
cat("MASE = ", knn_mase_fcd, "\n")
cat("SMAPE = ", knn_smape_fcd, "\n")
cat("RMSE = ", knn_rmse_fcd, "\n")
cat("R^2 = ", knn_r2_fcd)
```
