Mehrstufige Vorhersage mit dplyr und tun

Question

Die Do-Funktion in dplyr können Sie viele coole Modelle schnell und einfach machen, aber ich struggeling diese Modelle für gute rollende Prognosen zu verwenden.

# 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) 
        ) 

Mein Faktor variabel ist Hour, sind meine exogenen Variablen Wind und temp, und das, was ich will Price ist prognostizieren. Also, im Grunde habe ich 24 Modelle, mit denen ich rollende Vorhersagen machen könnte.

Jetzt enthält mein Datenrahmen 180 Tage. Ich würde gerne 100 Tage zurückgehen, und eine 1 Tag rollende Vorhersage machen und dann in der Lage sein, dies mit dem tatsächlichen Price zu vergleichen.

diese Brute-Force mache ich etwas würde wie folgt aussehen:

# 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 

Jetzt würde ich n ändern zu 99 * 24. Aber es wäre fantastisch, dies in einer Schleife haben oder anwenden, aber ich kann einfach‘ t herausfinden, wie es geht, und auch jede neue Prognose speichern.

Ich habe eine Schleife wie diese versucht, aber kein Glück noch:

# 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 
    } 

wirklich genial Bonus-Punkte für eine broom -ähnlichen Lösung :)

Answer 1

ich mit der Feststellung beginnen, werden feststellen, dass In Ihrer for-Schleife ist ein Fehler vorhanden. Anstelle von n*24*80 hast du wahrscheinlich (n+80)*24 gemeint. Der Zähler in Ihrer Schleife sollte auch von 0 bis 99 anstelle von 1 bis 100 gehen, wenn Sie auch die Vorhersage für den 81. Tag einbeziehen möchten.

Ich werde versuchen, eine elegante Lösung für Ihr Problem unten zu bieten. Zunächst definieren wir unsere Testdatenrahmen in genau der gleichen Art und Weise, die Sie in Ihrem Beitrag geleistet haben:

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) 
) 

Als nächstes definieren wir eine Funktion, die die Vorhersage für einen bestimmten Tag durchführt. Eingabeargumente sind der betrachtete Datenrahmen und die minimale Anzahl an Trainingstagen, die im Trainingset sein sollten (= 80 in diesem Beispiel). minTrainingDays+offSet+1 repräsentiert den tatsächlichen Tag, den wir vorhersagen. Beachten Sie, dass wir beginnen, für den Offset von 0 zu zählen.

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) 
    }) 
} 

Wir wollen Tage 81-180 voraussagen. Mit anderen Worten, wir benötigen mindestens 80 Tage in unserem Trainingssatz und möchten die Funktionsergebnisse für Offsets 0:99 berechnen. Dies kann mit einem einfachen lapply Aufruf erreicht werden. Wir beenden, indem Sie alle Ergebnisse zusammen in einem Datenrahmen verschmelzen:

# 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 Nach Ihrem Beitrag zu überprüfen und eine andere Antwort, die in der Zwischenzeit geschrieben wurde ich mit meiner Antwort zwei mögliche Probleme bemerkt. Zuerst wollten Sie ein rollendes Fenster von 80 Tagen Trainingsdaten. In meinem vorherigen Code werden jedoch alle verfügbaren Trainingsdaten verwendet, um das Modell anzupassen, anstatt nur 80 Tage zurückzugehen. Zweitens ist der Code nicht robust gegenüber Änderungen der Sommerzeit.

Diese beiden Probleme sind im folgenden Code behoben.Die Eingaben der Funktion sind jetzt auch intuitiver: Die Anzahl der Trainingstage und der tatsächliche vorhergesagte Tag können als Eingabeargumente verwendet werden. Beachten Sie, dass das Datenformat POSIXlt Dinge wie Sommerzeit, Schaltjahre usw. bei der Durchführung von Datenoperationen korrekt verarbeitet. Da die Daten in Ihrem Datenrahmen vom Typ POSIXct sind, müssen wir eine kleine Typumwandlung hin und her machen, um die Dinge richtig zu handhaben.

Neuer Code unten:

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) 

Ergebnisse wie folgt aussehen:

> 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

Man kann möglicherweise eine "rollende" data.frame mit dplyr erstellen, wie

folgt

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() 

Grundsätzlich , repliziert dieser Code die Beobachtungen nach Bedarf und weist jedem einen entsprechenden Day_Group_id zu 80-Tage-Stück. Dies ermöglicht es, group_by(Day_Group_id) zu verwenden, um das Modell separat für jeden 80-Tage-Chunk auszuführen.

Anschließend kann man es wie gewünscht verwenden. Kopieren Sie beispielsweise nur/Einfügen den arima Code von oben wie folgt:

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) 
    }) 

Bitte beachten Sie:

Die filter(n() >= 24*WINDOW_SIZE_DAYS - 1) verwendet wird hier statt filter(n() == 24*WINDOW_SIZE_DAYS) um volle 80-Tage-Fenster auszuwählen. Dies ist auf die Einstellung der Sommerzeit unter 2015-03-08 zurückzuführen. Die Stunde 2015-03-08 02:00:00 ist im Datensatz nicht vorhanden, da sie von 2015-03-08 01:00:00 direkt auf 2015-03-08 03:00:00 springt.