将滞后变量添加到lm模型中？

Question

我在一个时间序列上使用了lm，实际上效果非常好，并且速度超级超级快。

假设我的模型是：

> formula <- y ~ x

我在训练集上训练它：

> train <- data.frame( x = seq(1,3), y = c(2,1,4) )
> model <- lm( formula, train )

...而且我可以对新数据进行预测：

> test <- data.frame( x = seq(4,6) )
> test$y <- predict( model, newdata = test )
> test
  x        y
1 4 4.333333
2 5 5.333333
3 6 6.333333

这个功能表现非常好，速度也非常快。

我想在模型中添加滞后变量。我可以通过扩充我的原始训练集来实现这一点：

> train$y_1 <- c(0,train$y[1:nrow(train)-1])
> train
  x y y_1
1 1 2   0
2 2 1   2
3 3 4   1

更新公式：

formula <- y ~ x * y_1

...并且培训将会很顺利:

> model <- lm( formula, train )
> # no errors here

然而，问题在于没有办法使用“predict”，因为无法以批处理的方式填充测试集中的y_1。

现在，对于许多其他回归模型，有非常方便的方法来在公式中表达它们，例如poly(x,2)等，这些方法直接使用未经修改的训练和测试数据即可工作。

所以，我想知道是否有一种方法可以在公式中表达滞后变量，以便可以使用predict进行预测呢？理想情况下：

formula <- y ~ x * lag(y,-1)
model <- lm( formula, train )
test$y <- predict( model, newdata = test )

...不必增加（不确定是否是正确的词）训练和测试数据集，只需直接使用predict？

- Hugh Perkins

---

我认为R应该能够更加优雅地处理这个问题。 - Charlie

@Charlie，这个问题标记为“r”。你认为上面的代码是用什么语言编写的？ - Hugh Perkins

我知道它是用R编写的。我只是评论说，即使使用dyn包，我认为R在处理时间序列操作方面并不那么出色，并且我希望有一个可以更优雅地完成这项工作的包。例如，我认为Stata非常容易进行时间序列操作。 dyn包有助于回归，但是例如将滞后变量添加到数据框中需要一些技巧，如df$lagged <- c(NA, head(df$var, -1))。 - Charlie

啊，我明白了：这里的“should”是指“我希望它能做到”，而不是“我认为它应该做到”。 - Hugh Perkins

我认为你代码的最后一个块是有效的，前提是在覆盖之前test包含列y。 - user3226167

Answer 1

看看例如dynlm软件包，它可以给你提供滞后操作符。更一般地，Econometrics和Time Series的任务视图将有更多内容供您查看。

这是其示例的开头--一个月和十二个月的滞后：

R>      data("UKDriverDeaths", package = "datasets")
R>      uk <- log10(UKDriverDeaths)
R>      dfm <- dynlm(uk ~ L(uk, 1) + L(uk, 12))
R>      dfm

Time series regression with "ts" data:
Start = 1970(1), End = 1984(12)

Call:
dynlm(formula = uk ~ L(uk, 1) + L(uk, 12))

Coefficients:
(Intercept)     L(uk, 1)    L(uk, 12)  
      0.183        0.431        0.511  

R> 

- Dirk Eddelbuettel

---

如果您有时间，我该如何处理测试数据？例如，如果我在 train <- data.frame( y = head(UKDriverDeaths,96) ) 上进行训练，然后将我的测试数据设置为 test <- data.frame( y = rep(UKDriverDeaths[97],96) )，我会得到一条水平直线，即它使用了我的测试数据集中 y 的滞后值，而不是使用计算出的值。（编辑：使用 NA 也没有更好，即 test <- data.frame( y = c( UKDriverDeaths[97], rep(NA, 95) ) )：它只是给出了每个值的 NA）（编辑2：哦，也许可以使用 update？） - Hugh Perkins

似乎无法弄清如何进行预测。update(model,end=192) 似乎不起作用，model <- dynlm( y ~ L(y,1), end= 192) 也是如此。 - Hugh Perkins

我使用与dynlm库在语义上至少相关的dyn库尝试了一下，更新了问题。 - Hugh Perkins

Answer 2

在Dirk关于dynlm的建议下，我无法完全弄清楚如何进行预测，但是搜索这个问题让我发现了通过https://stats.stackexchange.com/questions/6758/1-step-ahead-predictions-with-dynlm-r-package使用dyn包的方法。

然后经过几个小时的尝试，我想出了以下函数来处理预测。在这个过程中有很多需要注意的事项，例如似乎不能rbind时间序列，而且predict的结果会受到start的影响以及其他一堆类似的问题，因此我认为这个答案相比仅仅命名一个包添加了很多价值，尽管我已经为Dirk的答案点赞。

因此，一个有效的解决方案是：

• 使用dyn包
• 使用以下方法进行预测

predictDyn方法：

# pass in training data, test data,
# it will step through one by one
# need to give dependent var name, so that it can make this into a timeseries
predictDyn <- function( model, train, test, dependentvarname ) {
    Ntrain <- nrow(train)
    Ntest <- nrow(test)
    # can't rbind ts's apparently, so convert to numeric first
    train[,dependentvarname] <- as.numeric(train[,dependentvarname])
    test[,dependentvarname] <- as.numeric(test[,dependentvarname])
    testtraindata <- rbind( train, test )
    testtraindata[,dependentvarname] <- ts( as.numeric( testtraindata[,dependentvarname] ) )
    for( i in 1:Ntest ) {
       result <- predict(model,newdata=testtraindata,subset=1:(Ntrain+i-1))
       testtraindata[Ntrain+i,dependentvarname] <- result[Ntrain + i + 1 - start(result)][1]
    }
    return( testtraindata[(Ntrain+1):(Ntrain + Ntest),] )
}

示例用法：

library("dyn")

# size of training and test data
N <- 6
predictN <- 10

# create training data, which we can get exact fit on, so we can check the results easily
traindata <- c(1,2)
for( i in 3:N ) { traindata[i] <- 0.5 + 1.3 * traindata[i-2] + 1.7 * traindata[i-1] }
train <- data.frame( y = ts( traindata ), foo = 1)

# create testing data, bunch of NAs
test <- data.frame( y = ts( rep(NA,predictN) ), foo = 1)

# fit a model
model <- dyn$lm( y ~ lag(y,-1) + lag(y,-2), train )
# look at the model, it's a perfect fit. Nice!
print(model)

test <- predictDyn( model, train, test, "y" )
print(test)

# nice plot
plot(test$y, type='l')

输出：

> model

Call:
lm(formula = dyn(y ~ lag(y, -1) + lag(y, -2)), data = train)

Coefficients:
(Intercept)   lag(y, -1)   lag(y, -2)  
        0.5          1.7          1.3  

> test
             y foo
7     143.2054   1
8     325.6810   1
9     740.3247   1
10   1682.4373   1
11   3823.0656   1
12   8686.8801   1
13  19738.1816   1
14  44848.3528   1
15 101902.3358   1
16 231537.3296   1

编辑：嗯，这个非常慢。即使我将subset中的数据限制为数据集的几行常数，每个预测大约需要24毫秒，对于我的任务来说，0.024*7*24*8*20*10/60/60=1.792小时 :-O

- Hugh Perkins

Answer 3

尝试使用ARIMA函数。AR参数是自回归，意味着滞后的y。xreg = 允许您添加其他X变量。您可以使用predict.ARIMA获得预测。

- PL84

Answer 4

有一个想法：

为什么不创建一个新的数据框？用所需的回归变量填充数据框。您可以为您想要的任何变量的所有滞后期设置列，例如L1、L2、...、Lp，然后您可以像进行交叉类型回归一样使用您的函数。

因为您不必每次调用拟合和预测函数时都操作数据，而是只需转换一次数据，因此速度会更快。我知道Eviews和Stata提供滞后运算符。这确实很方便。但是，如果您不需要像“lm”函数计算的所有内容，则效率也很低。如果您需要执行数十万次迭代，并且只需要预测或预测以及诸如BIC或AIC之类的信息准则的值，则通过避免进行不使用的计算，您可以通过编写OLS估计器函数来击败“lm”在速度上。

- Stéphane

---

你知道在行中添加滞后值作为列的实用方法吗？我认为我们需要一个所需滞后阶数的窗口，并将其通过日期进行移动以获取滞后值和相应的输出值。 - ibilgen