Previsione a più passaggi con dplyr e do

Question

La funzione di do in dplyr consente di creare molti modelli interessanti in modo facile e veloce, ma mi sto sforzando di utilizzare questi modelli per le buone previsioni rolling.

# Data illustration 

require(dplyr) 
require(forecast) 

df <- data.frame(
    Date = seq.POSIXt(from = as.POSIXct("2015-01-01 00:00:00"), 
        to = as.POSIXct("2015-06-30 00:00:00"), by = "hour")) 

    df <- df %>% mutate(Hour = as.numeric(format(Date, "%H")) + 1, 
         Wind = runif(4320, min = 1, max = 5000), 
         Temp = runif(4320, min = - 20, max = 25), 
         Price = runif(4320, min = -15, max = 45) 
        ) 

mio variabile fattore è Hour, le mie variabili esogene sono Wind e temp, e la cosa che voglio previsione è Price. Quindi, in sostanza, ho 24 modelli che mi piacerebbe essere in grado di fare previsioni a rotazione con.

Ora il mio frame di dati contiene 180 giorni. Vorrei tornare indietro di 100 giorni e fare una previsione a rotazione di 1 giorno per poterlo confrontare con l'attuale Price.

In questo modo la forza bruta sarebbe simile a questa:

# First I fit the data frame to be exactly the right length 
# 100 days to start with (2015-03-21 or so), then 99, then 98.., etc. 
n <- 100 * 24 

# Make the price <- NA so I can replace it with a forecast 
df$Price[(nrow(df) - n): (nrow(df) - n + 24)] <- NA 

# Now I make df just 81 days long, the estimation period + the first forecast 
df <- df[1 : (nrow(df) - n + 24), ] 

# The actual do & fit, later termed fx(df) 

result <- df %>% group_by(Hour) %>% do ({ 
    historical <- .[!is.na(.$Price), ] 
    forecasted <- .[is.na(.$Price), c("Date", "Hour", "Wind", "Temp")] 
    fit <- Arima(historical$Price, xreg = historical[, 3:4], order = c(1, 1, 0)) 
    data.frame(forecasted[], 
      Price = forecast.Arima(fit, xreg = forecasted[3:4])$mean) 
}) 

result 

Ora vorrei cambiare n-99 * 24. Ma sarebbe fantastico avere questo in un ciclo o applicare, ma ho semplicemente puo' t capire come farlo, e anche salvare ogni nuova previsione.

Ho provato un circuito come questo, ma ancora senza fortuna:

# 100 days ago, forecast that day, then the next, etc. 
for (n in 1:100) { 
    nx <- n * 24 * 80   # Because I want to start after 80 days 
    df[nx:(nx + 23), 5] <- NA # Set prices to NA so I can forecast them 
    fx(df) # do the function 
    df.results[n] <- # Write the results into a vector/data frame to save them 
    # and now rinse and repeat for n + 1 
    } 

veramente impressionante bonus-punti per una soluzione :) broom -come

Answer 1

Comincerò osservando che c'è un errore presente nel tuo ciclo for. Invece di n*24*80 probabilmente intendevi (n+80)*24. Il contatore del ciclo dovrebbe anche andare da 0 a 99 anziché da 1 a 100 se si desidera includere anche la previsione per l'81 ° giorno.

Proverò a fornire una soluzione elegante per il tuo problema di seguito. In primo luogo, definiamo il nostro dataframe test nello stesso modo esatto in cui hai fatto nel tuo post:

set.seed(2) 
df <- data.frame(
Date = seq.POSIXt(from = as.POSIXct("2015-01-01 00:00:00"), 
        to = as.POSIXct("2015-06-30 00:00:00"), by = "hour")) 
df <- df %>% mutate(Hour = as.numeric(format(Date, "%H")) + 1, 
        Wind = runif(4320, min = 1, max = 5000), 
        Temp = runif(4320, min = - 20, max = 25), 
        Price = runif(4320, min = -15, max = 45) 
) 

successivo, definiamo una funzione che esegue la previsione per un giorno particolare. Gli argomenti di input sono il dataframe considerato e il numero minimo di giorni di formazione che dovrebbero essere nel set di allenamento (= 80 in questo esempio). minTrainingDays+offSet+1 rappresenta il giorno reale che stiamo prevedendo. Si noti che iniziamo a contare da 0 per l'offset.

forecastOneDay <- function(theData,minTrainingDays,offset) 
{ 
    nrTrainingRows <- (minTrainingDays+offset)*24 

    theForecast <- theData %>% 
    filter(min_rank(Date) <= nrTrainingRows+24) %>% # Drop future data that we don't need 
    group_by(Hour) %>% 
    do ({ 
     trainingData <- head(.,-1) # For each group, drop the last entry from the dataframe 
     forecastData <- tail(.,1) %>% select(Date,Hour,Wind,Temp) # For each group, predict the last entry 
     fit <- Arima(trainingData$Price, xreg=trainingData[,3:4], order=c(1,1,0)) 
     data.frame(forecastData, realPrice = tail(.,1)$Price, predictedPrice = forecast.Arima(fit,xreg=forecastData[3:4])$mean) 
    }) 
} 

Vogliamo prevedere i giorni 81-180. In altre parole, abbiamo bisogno di un minimo di 80 giorni nel nostro set di allenamento e vogliamo calcolare i risultati delle funzioni per gli offset 0:99. Questo può essere fatto con una semplice chiamata lapply. Finiamo fondendo tutti i risultati insieme in una dataframe:

# Perform one day forecasts for days 81-180 
resultList <- lapply(0:99, function(x) forecastOneDay(df,80,x)) 
# Merge all the results 
mergedForecasts <- do.call("rbind",resultList) 

EDIT Dopo aver esaminato il tuo post e un'altra risposta che è stato pubblicato nel frattempo ho notato due potenziali problemi con la mia risposta. In primo luogo, si desiderava una a rotazione della finestra di 80 giorni di dati di allenamento. Tuttavia, nel mio codice precedente tutti i dati di allenamento disponibili vengono utilizzati per adattarsi al modello invece di tornare indietro di soli 80 giorni. In secondo luogo, il codice non è robusto alle modifiche dell'ora legale.

Questi due problemi sono stati risolti nel seguente codice.Gli input della funzione sono anche più intuitivi ora: il numero di giorni di formazione e il giorno previsto effettivo possono essere utilizzati come parametri di input. Si noti che il formato dati POSIXlt gestisce le cose come l'ora legale, gli anni bisestili ecc. Correttamente quando si eseguono operazioni sulle date. Poiché le date nel tuo dataframe sono di tipo POSIXct, dobbiamo eseguire una piccola conversione di tipo avanti e indietro per gestire le cose correttamente.

Nuovo codice qui sotto:

forecastOneDay <- function(theData,nrTrainingDays,predictDay) # predictDay should be greater than nrTrainingDays 
{ 
    initialDate <- as.POSIXlt(theData$Date[1]); # First day (midnight hour) 
    startDate <- initialDate # Beginning of training interval 
    endDate <- initialDate # End of test interval 

    startDate$mday <- initialDate$mday + (predictDay-nrTrainingDays-1) # Go back 80 days from predictday 
    endDate$mday <- startDate$mday + (nrTrainingDays+1) # +1 to include prediction day 

    theForecast <- theData %>% 
    filter(Date >= as.POSIXct(startDate),Date < as.POSIXct(endDate)) %>% 
    group_by(Hour) %>% 
    do ({ 
     trainingData <- head(.,-1) # For each group, drop the last entry from the dataframe 
     forecastData <- tail(.,1) %>% select(Date,Hour,Wind,Temp) # For each group, predict the last entry 
     fit <- Arima(trainingData$Price, xreg=trainingData[,3:4], order=c(1,1,0)) 
     data.frame(forecastData, realPrice = tail(.,1)$Price, predictedPrice = forecast.Arima(fit,xreg=forecastData[3:4])$mean) 
    }) 
} 

# Perform one day forecasts for days 81-180 
resultList <- lapply(81:180, function(x) forecastOneDay(df,80,x)) 
# Merge all the results 
mergedForecasts <- do.call("rbind",resultList) 

risultati simile a questa:

> head(mergedForecasts) 
Source: local data frame [6 x 6] 
Groups: Hour 

       Date Hour  Wind  Temp realPrice predictedPrice 
1 2015-03-22 00:00:00 1 1691.589 -8.722152 -11.207139  5.918541 
2 2015-03-22 01:00:00 2 1790.928 18.098358 3.902686  37.885532 
3 2015-03-22 02:00:00 3 1457.195 10.166422 22.193270  34.984164 
4 2015-03-22 03:00:00 4 1414.502 4.993783 6.370435  12.037642 
5 2015-03-22 04:00:00 5 3020.755 9.540715 25.440357  -1.030102 
6 2015-03-22 05:00:00 6 4102.651 2.446729 33.528199  39.607848 
> tail(mergedForecasts) 
Source: local data frame [6 x 6] 
Groups: Hour 

       Date Hour  Wind  Temp realPrice predictedPrice 
1 2015-06-29 18:00:00 19 1521.9609 13.6414797 12.884175  -6.7789109 
2 2015-06-29 19:00:00 20 555.1534 3.4758159 37.958768  -5.1193514 
3 2015-06-29 20:00:00 21 4337.6605 4.7242352 -9.244882  33.6817379 
4 2015-06-29 21:00:00 22 3140.1531 0.8127839 15.825230  -0.4625457 
5 2015-06-29 22:00:00 23 1389.0330 20.4667234 -14.802268  15.6755880 
6 2015-06-29 23:00:00 24 763.0704 9.1646139 23.407525  3.8214642 

Answer 2

Uno potenzialmente in grado di creare un data.frame "rolling" con dplyr come segue

library(dplyr) 
library(lubridate) 

WINDOW_SIZE_DAYS <- 80 

df2 <- df %>% 
    mutate(Day = yday(Date)) %>% 
    replicate(n = WINDOW_SIZE_DAYS, simplify = FALSE) %>% 
    bind_rows %>% 
    group_by(Date) %>% 
    mutate(Replica_Num = 1:n()) %>% 
    mutate(Day_Group_id = Day + Replica_Num - 1) %>% 
    ungroup() %>% 
    group_by(Day_Group_id) %>% 
    filter(n() >= 24*WINDOW_SIZE_DAYS - 1) %>% 
    select(-Replica_Num) %>% 
    arrange(Date) %>% 
    ungroup() 

Fondamentalmente , questo codice replica le osservazioni secondo necessità e assegna un corrispondente Day_Group_id a ciascuna Pezzo di 80 giorni. Ciò che consente è di utilizzare group_by(Day_Group_id) per eseguire il modello separatamente su ogni blocco di 80 giorni.

Successivamente, è possibile utilizzarlo come desiderato. Ad esempio, basta copiare/incollare il codice Arima dall'alto come segue:

df3 <- df2 %>% 
    group_by(Day_Group_id, Hour) %>% 
    arrange(Date) %>% 
    do ({ 
    trainingData <- head(.,-1) # For each group, drop the last entry from the dataframe 
    forecastData <- tail(.,1) %>% select(Date,Hour,Wind,Temp) # For each group, predict the last entry 
    fit <- Arima(trainingData$Price, xreg=trainingData[,3:4], order=c(1,1,0)) 
    data.frame(forecastData, realPrice = tail(.,1)$Price, predictedPrice = forecast.Arima(fit,xreg=forecastData[3:4])$mean) 
    }) 

Attenzione: per

Il filter(n() >= 24*WINDOW_SIZE_DAYS - 1) è usato qui al posto di filter(n() == 24*WINDOW_SIZE_DAYS) al fine di selezionare le finestre piene di 80 giorni. Ciò è dovuto alla regolazione dell'ora legale su 2015-03-08. L'ora 2015-03-08 02:00:00 non esiste nel set di dati mentre salta da 2015-03-08 01:00:00 direttamente allo 2015-03-08 03:00:00.