r - 使用 ggplot2 和 funggcast 函数进行预测

Question

在 this在网站上，Davenport 先生发布了一个函数，可以使用 ggplot2 绘制有马预测。以任意数据集为例，他发表了 here .我可以跟随他的例子而没有任何错误信息。

现在，当我使用我的数据时，我会以警告结束:

1: In window.default(x, ...) : 'end' value not changed
2: In window.default(x, ...) : 'end' value not changed

我知道当我调用这个命令时会发生这种情况 pd <- funggcast(yt, yfor)由于我在数据中指出的数据存在问题 end = c(2013) .但我不知道如何解决这个问题。

这是我使用的代码:

library(ggplot2)
library(zoo)
library(forecast)

myts <- ts(rnorm(55), start = c(1960), end = c(2013), freq = 1)
funggcast <- function(dn, fcast){

en <- max(time(fcast$mean)) # Extract the max date used in the forecast

# Extract Source and Training Data
ds <- as.data.frame(window(dn, end = en))
names(ds) <- 'observed'
ds$date <- as.Date(time(window(dn, end = en)))

# Extract the Fitted Values (need to figure out how to grab confidence intervals)
dfit <- as.data.frame(fcast$fitted)
dfit$date <- as.Date(time(fcast$fitted))
names(dfit)[1] <- 'fitted'

ds <- merge(ds, dfit, all.x = T) # Merge fitted values with source and training data

# Extract the Forecast values and confidence intervals
dfcastn <- as.data.frame(fcast)
dfcastn$date <- as.Date(as.yearmon(row.names(dfcastn)))
names(dfcastn) <- c('forecast','lo80','hi80','lo95','hi95','date')

pd <- merge(ds, dfcastn,all.x = T) # final data.frame for use in ggplot
return(pd)

} 

yt <- window(myts, end = c(2013)) # extract training data until last year
yfit <- auto.arima(myts) # fit arima model
yfor <- forecast(yfit) # forecast
pd <- funggcast(yt, yfor) # extract the data for ggplot using function funggcast()

ggplot(data = pd, aes(x = date,y = observed)) + geom_line(color = "red") + geom_line(aes(y = fitted), color = "blue") + geom_line(aes(y = forecast)) + geom_ribbon(aes(ymin = lo95, ymax = hi95), alpha = .25) + scale_x_date(name = "Time in Decades") + scale_y_continuous(name = "GDP per capita (current US$)") + theme(axis.text.x = element_text(size = 10), legend.justification=c(0,1), legend.position=c(0,1)) + ggtitle("Arima(0,1,1) Fit and Forecast of GDP per capita for Brazil (1960-2013)") + scale_color_manual(values = c("Blue", "Red"), breaks = c("Fitted", "Data", "Forecast"))

编辑:我发现了另一个博客 here具有与 forecast 一起使用的功能和 ggplot2但如果我能找到我的错误，我想使用上面的方法。任何人？

编辑 2:
如果我用我的数据运行你更新的代码 here ，然后我把图表放在下面。请注意，我没有更改 end = c(2023)为 mtys ，否则不会将预测值与拟合值合并。

myts <- ts(WDI_gdp_capita$Brazil, start = c(1960), end = c(2023), freq = 1)

funggcast <- function(dn, fcast){

  en <- max(time(fcast$mean)) # Extract the max date used in the forecast

  # Extract Source and Training Data
  ds <- as.data.frame(window(dn, end = en))
  names(ds) <- 'observed'
  ds$date <- as.Date(time(window(dn, end = en)))

  # Extract the Fitted Values (need to figure out how to grab confidence intervals)
  dfit <- as.data.frame(fcast$fitted)
  dfit$date <- as.Date(time(fcast$fitted))
  names(dfit)[1] <- 'fitted'

  ds <- merge(ds, dfit, all = T) # Merge fitted values with source and training data

  # Extract the Forecast values and confidence intervals
  dfcastn <- as.data.frame(fcast)
  dfcastn$date <- as.Date(paste(row.names(dfcastn),"01","01",sep="-"))
  names(dfcastn) <- c('forecast','lo80','hi80','lo95','hi95','date')

  pd <- merge(ds, dfcastn,all.x = T) # final data.frame for use in ggplot
  return(pd)

} # ggplot function by Frank Davenport

yt <- window(myts, end = c(2013)) # extract training data until last year
yfit <- auto.arima(yt) # fit arima model
yfor <- forecast(yfit) # forecast
pd <- funggcast(myts, yfor) # extract the data for ggplot using function funggcast()

ggplot(data = pd, aes(x = date, y = observed)) + geom_line(color = "red") + geom_line(aes(y = fitted), color = "blue") + geom_line(aes(y = forecast)) + geom_ribbon(aes(ymin = lo95, ymax = hi95), alpha = .25) + scale_x_date(name = "Time in Decades") + scale_y_continuous(name = "GDP per capita (current US$)") + theme(axis.text.x = element_text(size = 10), legend.justification=c(0,1), legend.position=c(0,1)) + ggtitle("Arima(0,1,1) Fit and Forecast of GDP per capita for Brazil (1960-2013)") + scale_color_manual(values = c("Blue", "Red"), breaks = c("Fitted", "Data", "Forecast")) + ggsave((filename = "gdp_forecast_ggplot.pdf"), width=330, height=180, units=c("mm"), dpi = 300, limitsize = TRUE)

我得到的几乎完美的图表:


另一个问题:如何在此图中获得图例？

如果我设置 end = c(2013)为 myts ，我得到与开始时相同的图形:

Answer 1

达文波特先生的分析与您试图制作的情节之间有几点不同。
第一个是他将 arima 预测与一些观察到的数据进行比较，这就是他在整个时间序列的一部分(训练集)上训练模型的原因。
为此，您应该使初始时间序列更长:

myts <- ts(rnorm(55), start = c(1960), end = c(2023), freq = 1)

然后在脚本的末尾，您可以选择截至 2013 年的培训:

yt <- window(myts, end = c(2013)) # extract training data until last year

模型应该在训练集上训练，而不是整个时间序列，因此您应该将 yfit 行更改为:

yfit <- auto.arima(yt) # fit arima model

并使用整个时间序列调用 funggcast 函数，因为它需要观察和拟合的数据:

pd <- funggcast(myts, yfor)

最后，他使用具有月份和年份的日期，因此在他的 funggcast 中功能，更改这一行:

dfcastn$date <- as.Date(as.yearmon(row.names(dfcastn)))

到:

dfcastn$date <- as.Date(paste(row.names(dfcastn),"01","01",sep="-"))

这是因为模型预测的值需要更改为日期，例如必须将 2014 更改为 2014-01-01，以便与观察到的数据合并。

全部更改后，代码如下所示:

library(ggplot2)
library(zoo)
library(forecast)

myts <- ts(rnorm(55), start = c(1960), end = c(2013), freq = 1)
funggcast <- function(dn, fcast){

        en <- max(time(fcast$mean)) # Extract the max date used in the forecast

        # Extract Source and Training Data
        ds <- as.data.frame(window(dn, end = en))
        names(ds) <- 'observed'
        ds$date <- as.Date(time(window(dn, end = en)))

        # Extract the Fitted Values (need to figure out how to grab confidence intervals)
        dfit <- as.data.frame(fcast$fitted)
        dfit$date <- as.Date(time(fcast$fitted))
        names(dfit)[1] <- 'fitted'

        ds <- merge(ds, dfit, all.x = T) # Merge fitted values with source and training data

        # Extract the Forecast values and confidence intervals
        dfcastn <- as.data.frame(fcast)
        dfcastn$date <- as.Date(paste(row.names(dfcastn),"01","01",sep="-"))
        names(dfcastn) <- c('forecast','lo80','hi80','lo95','hi95','date')

        pd <- merge(ds, dfcastn,all= T) # final data.frame for use in ggplot
        return(pd)

} 

yt <- window(myts, end = c(2013)) # extract training data until last year
yfit <- auto.arima(yt) # fit arima model
yfor <- forecast(yfit) # forecast
pd <- funggcast(myts, yfor) # extract the data for ggplot using function funggcast()

plotData<-ggplot(data = pd, aes(x = date, y = observed)) + geom_line(aes(color = "1")) +
        geom_line(aes(y = fitted,color="2")) + 
        geom_line(aes(y = forecast,color="3")) +
        scale_colour_manual(values=c("red", "blue","black"),labels = c("Observed", "Fitted", "Forecasted"),name="Data")+
        geom_ribbon(aes(ymin = lo95, ymax = hi95), alpha = .25)+
        scale_x_date(name = "Time in Decades") +
        scale_y_continuous(name = "GDP per capita (current US$)")+
        theme(axis.text.x = element_text(size = 10)) + 
        ggtitle("Arima(0,1,1) Fit and Forecast of GDP per capita for Brazil (1960-2013)")

plotData

你会得到一个看起来像这样的图，完全随机的时间序列的拟合非常糟糕。另外 ggplot 会输出一些错误，因为预测线在 2013 年之前没有数据，并且拟合数据在 2013 年之后没有继续。(我运行了几次，根据初始的随机时间序列，模型可能只是在任何地方预测 0)



编辑:更改了 pd分配线也是如此，以防在 2013 年之后没有观察到的数据

Edit2:我在代码末尾更改了 ggplot 函数以确保图例显示