【发布时间】：2021-09-14 01:06:03
【问题描述】：
我正在管理包含 348 次观察的纯月度时间序列数据。
这是可重现的样本数据：
library(dplyr)
set.seed(123)
Y <- cumsum(rnorm(48))
date <- as.Date(c("2012-01-01", "2012-02-01", "2012-03-01", "2012-04-01",
"2012-05-01","2012-06-01", "2012-07-01", "2012-08-01",
"2012-09-01","2012-10-01","2012-11-01", "2012-12-01",
"2013-01-01", "2013-02-01","2013-03-01", "2013-04-01",
"2013-05-01","2013-06-01", "2013-07-01", "2013-08-01",
"2013-09-01","2013-10-01","2013-11-01", "2013-12-01",
"2014-01-01", "2014-02-01","2014-03-01", "2014-04-01",
"2014-05-01","2014-06-01", "2014-07-01", "2014-08-01",
"2014-09-01","2014-10-01","2014-11-01", "2014-12-01",
"2015-01-01", "2015-02-01", "2015-03-01", "2015-04-01",
"2015-05-01","2015-06-01", "2015-07-01", "2015-08-01",
"2015-09-01","2015-10-01","2015-11-01", "2015-12-01"))
data<-data.frame(date,Y)
我正在复制一篇论文并计算“冲击”。引用的程序如下：
“每个系列中的冲击均由 AR(2) 模型在 10 个月的滚动窗口内计算得出，该滚动窗口在第 n 个月结束。第 n+1 个月的冲击表示为 dYn+1 是实际值之间的差异使用过去 10 个月估计的斜率系数计算系列及其预测值。因此，我们的方法具有前瞻性，提供了样本外预测误差。"
原作者提出的带有趋势项的AR(2)模型如下：
Yt = a0 + a1*Yt-1 + a2*Yt-2 + a3*Tt+ residualt
where Tt is the serial number of the observation, to account for a time trend in these series.
如果我的目标是使用前 10 个 obs 计算第 11 个中的平均值，我可以简单地调用以下代码：
Mean= slide_index_dbl(Y, date, mean, .before = months(10), .after = months(-1), .complete = T)
但是，在这种情况下，目标是使用所有之前的 10 个 obs 运行 AR 模型，并使用此估计模型来预测第 11 个，最终输出是第 11 个中的实际值减去预测值。简单地说，我需要构造一个函数来实现这个目标，而不是在前面的例子中使用“mean”函数。
一旦我完成了这个函数（我们称之为 AR_2），我们就可以在滑块内部调用它。
library(slider)
data1<-data%>%
mutate(Shock= slide_index_dbl(Y, date, AR_2, .before = months(10), .after = months(-1), .complete = T))
Date Y N Predict Shock
2012-01-01 0.15 1 0.2 -0.005
2012-02-01 0.4 2 0.33 0.07
2012-03-01 0.39 3 0.44 -0.05
...
2012-10-01 1.85 10 2.1 -0.25
2012-11-01 1.7 11 1.5 0.2
2012-12-01 3.46 12 4.1 -0.65
让我使用上面我编写的示例数据来说明我的问题。最终的输出是 Shock，它是 Y（实际数据）和 Predicted Y 之间的差异。话虽如此，问题是如何得到 Predicted Y。在我们预测 Y 之前，我们需要先训练一个 AR 模型，通过输入先前的10 个月的数据。一旦你得到这个模型，你就可以用这个模型预测第 11 个 obs，这将是 Predicted Y。最后的尝试是计算它与实际 Y 之间的差异，这称为冲击。
数值示例如下：为了获得 2012-11-01 的输出 0.074，我需要使用从 2012-01-01 到 2012-10- 的所有前 10 个月数据训练一个 AR(2) 模型01，即 Y 的 0.15 到 1.85 和 N 的 1 到 10。一旦我训练了这个模型，我用它来预测 2012-11-01 的新值（0.2）。最终输出“Shock”是 2012-11-01 实际 obs 与预测 obs 之间的差异 (0.15-0.2=0.05)。
同样，为了获得 2012-12-01 的输出 0.08，我需要使用从 2012-02-01 到 2012-11-01 的所有前 10 个月数据训练一个 AR(2) 模型，即来自Y 为 0.4 到 1.7，N 为 2 到 11。一旦我训练了这个模型，我就会用它来预测 2012 年 12 月 1 日的新值 (4.1)。最终输出“Shock”是 2012-12-01 实际 obs 与预测 obs 之间的差异 (3.46-4.1=-0.65)。
我不知道如何编写这样的函数（AR_2）并在 slide_index_dbl 中调用它。请记住 AR_2 中有一个趋势项 a3*Tt，我不知道如何建模。
以下是我尝试过的。我使用 lm 来实现 AR moel 而不是 arima。这是因为我不知道如何使用 arima 来预测新值。如果有人熟悉 arima 或 Arima 功能，请使用它。
Intercept_extract_lm<-function(x){
N<-rep(1:10)
model<-lm(x~ lag(x,1)+ lag(x, 2)+N)
coef(model)["(Intercept)"]
}
Log_1_extract_lm<-function(x){
N<-rep(1:10)
model<-lm(x~ lag(x,1)+ lag(x, 2)+N)
coef(model)["Lag_1"]
}
Log_2_extract_lm<-function(x){
N<-rep(1:10)
model<-lm(x~ lag(x,1)+ lag(x, 2)+N)
coef(model)["Lag_2"]
}
drift_extract_lm<-function(x){
N<-rep(1:10)
model<-lm(x~ lag(x,1)+ lag(x, 2)+N)
coef(model)["N"]
}
data1<-data%>%
mutate(Lag_1=lag(Y,1),Lag_2=lag(Y,2),N=1:n(),
a1=slide_index_dbl(Y, Date, Log_1_extract_lm, .before = months(10), .after = months(-1), .complete = T),
a2=slide_index_dbl(Y, Log_2_extract_lm, .before = months(10), .after = months(-1), .complete = T),
drift=slide_index_dbl(Y, Date, drift_extract_lm, .before = months(10), .after = months(-1), .complete = T),
Intercept = slide_index_dbl(Y, Date, Intercept_extract_lm, .before = months(10), .after = months(-1), .complete = T),
Predict=Intercept+a1*Lag_1+a2*Lag_2+drift,
Shock=Y-Predict)
我知道我的代码中最大的问题是它只能在滑块中接受一个输入参数，但我在定义的所有自定义函数中使用了四个（Y 和 lag_1 和 lag_2 和 N）。与前面的例子不同，我们只计算一个变量的滚动平均值，在这种情况下，为了运行这样的回归，我们需要在每个滚动窗口中有四个变量，但滑块只有一个变量输入。即使我们编辑“提取”函数将四个变量减少为两个（lag_1 和 lag_2 可以写为 lag(Y,1) 和 lag(Y,2)，我们仍然有 Y 和 N 两个变量输入）
关于趋势项，据我了解，例如2012-11-01，为了运行回归，我需要前10分钟的N，即1,2,3....10，即是obs的序列号。对于 2012 年 12 月 1 日，它还应该包括 1,2,3...10 中的 N，但是，根据我的代码，它是 2,3,4...11。但是添加常数 1 (2,3...11 vs 1,2,...10) 不会影响回归系数，对吧？老实说，我不确定这种情况下的趋势术语，如果您了解原始作者，请随时更改它。
【问题讨论】：
-
嘿，阿尼尔。我已经编辑了这个问题。这是我在另一个问题中提到的问题。最初的程序是滚动 60 个月，但考虑到我的样本有限，我改为 10 个月。每个月最终输出的长度只有1。我只是不知道如何在slider中编写这样的AR_2函数。
-
阿尼尔，我回来看看你对这个问题有什么想法吗？我在上一个问题中看到了您的解释，您使用 mean(arima(x,c（1,0，0）0$residuals) 来提取长度为 1 的值。但是，当我返回程序时，似乎这不是原作者想要的。
-
阿尼尔，请看我编辑的问题。这个问题可能需要一些时间序列分析背景，因为我将使用 ARIMA（AR 模型）建模。让我知道你的困惑。
-
阿尼尔，我已经编辑了最后几段。如果您有在 R 中模拟 ARIMA 模型的经验，那会更容易。困难在于使用由前 10 个 obs 提供的估计 ARIMA 模型来预测新值。当我们要计算前一个obs的平均值时，其余部分与上一个问题相同。
-
阿尼尔，我卸载了我尝试过的东西。我明白我的错误意味着什么，我只是不知道如何解决它。因为在滑块中，变量是 Y，所以我只使用前 10 个 Y，但是，在回归中，我还需要使用前 10 个 Lag_1 和 Lag_2。
标签： r regression arima rolling-computation forecast