R 中的滚动 AR 回归
Rolling AR regressions in R
我正在管理包含 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) 模型计算的，滚动 window 10 个月，结束于 n.The 月 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.
如果我的 objective 是使用前 10 个 obs 计算 11st 的平均值，我可以简单地调用以下内容：
Mean= slide_index_dbl(Y, date, mean, .before = months(10), .after = months(-1), .complete = T)
然而，在这种情况下，objective是运行一个AR模型，使用所有之前的10个obs并使用这个估计模型来预测第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(actual data)和Predicted Y的差值。话虽如此，问题是如何得到Predicted Y。在我们预测Y之前，我们需要先训练一个AR模型，通过feed previous 10个月的数据。一旦你得到这个模型，你就可以用这个模型预测11st obs，那就是Predicted Y。最后的尝试是计算它与实际Y之间的差异，这称为shock。
数值示例如 follows:in 为了获得 2012-11-01 的输出 0.074，我需要使用从 2012-01-01 到2012-10-01，即 Y 从 0.15 到 1.85，N 从 1 到 10。训练完这个模型后，我用它来预测 2012-11-01 的新值 (0.2)。最终输出“Shock”是 2012-11-01 实际 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 与预测值 (3.46-4.1=-0.65) 之间的差异。
我不知道如何编写这样的函数(AR_2)并在slide_index_dbl中调用它。请记住 AR_2 中有一个趋势项 a3*Tt，我不确定如何对其建模。
以下是我试过的。我使用 lm 来实现 AR 模型而不是 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）。不像前面的例子，我们只是计算一个单一变量的滚动平均值，在这种情况下，为了运行这样的回归，我们需要在每个滚动window中有四个变量，但是滑块只有一个变量输入.即使我们编辑“提取”函数将四个变量减少为两个（lag_1 和 lag_2 可以写成 lag(Y,1) 和 lag(Y,2)，我们仍然有 Y 和 N两个变量输入)
关于趋势项，据我了解，例如2012-11-01，为了运行回归，我需要前10分钟N，即1,2,3。 ...10，即obs的序号。对于 2012-12-01，它还应该包括来自 1,2,3...10 的 N，但是，根据我的代码，它是 2,3,4...11。但是添加常数 1（2,3...11 对 1,2,...10）不会影响回归系数，对吧？老实说，我不太确定这种情况下的趋势术语，如果您了解原作者，请随时更改它。
你可以这样做。您的滞后 2 [A(2) 模型] 问题将由 order = c(2,0,0) 解决，您不必明确地这样做。 arima() 函数与 predict().
结合使用
library(tidyverse)
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)
library(lubridate, warn.conflicts = F)
library(slider)
data %>%
mutate(pred_roll_10 = slide_index_dbl(Y, .i = date,
~ predict(arima(.x, c(2,0,0), method = 'ML'),1)$pred[1],
.before = months(10),
.after = months(-1),
.complete = T),
shock = Y - pred_roll_10)
#> Warning in log(s2): NaNs produced
#> Warning in log(s2): NaNs produced
#> Warning in log(s2): NaNs produced
#> Warning in log(s2): NaNs produced
#> date Y pred_roll_10 shock
#> 1 2012-01-01 -0.5604756 NA NA
#> 2 2012-02-01 -0.7906531 NA NA
#> 3 2012-03-01 0.7680552 NA NA
#> 4 2012-04-01 0.8385636 NA NA
#> 5 2012-05-01 0.9678513 NA NA
#> 6 2012-06-01 2.6829163 NA NA
#> 7 2012-07-01 3.1438325 NA NA
#> 8 2012-08-01 1.8787713 NA NA
#> 9 2012-09-01 1.1919184 NA NA
#> 10 2012-10-01 0.7462564 NA NA
#> 11 2012-11-01 1.9703382 0.6951079 1.275230339
#> 12 2012-12-01 2.3301521 2.0524671 0.277684986
#> 13 2013-01-01 2.7309235 1.9333295 0.797594017
#> 14 2013-02-01 2.8416062 2.0486834 0.792922863
#> 15 2013-03-01 2.2857651 1.9959595 0.289805588
#> 16 2013-04-01 4.0726782 1.7514948 2.321183420
#> 17 2013-05-01 4.5705287 3.6469980 0.923530738
#> 18 2013-06-01 2.6039116 4.0585260 -1.454614411
#> 19 2013-07-01 3.3052675 1.6480334 1.657234048
#> 20 2013-08-01 2.8324760 2.9543635 -0.121887466
#> 21 2013-09-01 1.7646523 2.8739943 -1.109341911
#> 22 2013-10-01 1.5466774 2.7845341 -1.237856702
#> 23 2013-11-01 0.5206730 2.5894221 -2.068749114
#> 24 2013-12-01 -0.2082182 1.3652172 -1.573435497
#> 25 2014-01-01 -0.8332575 0.3654956 -1.198753080
#> 26 2014-02-01 -2.5199508 -0.5625776 -1.957373251
#> 27 2014-03-01 -1.6821638 -2.6388755 0.956711703
#> 28 2014-04-01 -1.5287907 -0.7131115 -0.815679154
#> 29 2014-05-01 -2.6669276 -1.2140023 -1.452925253
#> 30 2014-06-01 -1.4131127 -2.6421098 1.228997135
#> 31 2014-07-01 -0.9866485 -1.0140132 0.027364706
#> 32 2014-08-01 -1.2817199 -0.8209904 -0.460729511
#> 33 2014-09-01 -0.3865943 -1.2655968 0.879002504
#> 34 2014-10-01 0.4915392 -1.1977841 1.689323273
#> 35 2014-11-01 1.3131203 -0.4463373 1.759457622
#> 36 2014-12-01 2.0017605 0.9931931 1.008567426
#> 37 2015-01-01 2.5556782 1.7819117 0.773766509
#> 38 2015-02-01 2.4937665 2.3699293 0.123837184
#> 39 2015-03-01 2.1878038 2.1804620 0.007341787
#> 40 2015-04-01 1.8073328 1.5452418 0.262091026
#> 41 2015-05-01 1.1126258 1.2644697 -0.151843849
#> 42 2015-06-01 0.9047086 0.2974574 0.607251157
#> 43 2015-07-01 -0.3606878 0.7761950 -1.136882824
#> 44 2015-08-01 1.8082682 -1.1719233 2.980191448
#> 45 2015-09-01 3.0162302 1.9584189 1.057811310
#> 46 2015-10-01 1.8931216 2.5356767 -0.642555104
#> 47 2015-11-01 1.4902368 1.2115509 0.278685838
#> 48 2015-12-01 1.0235814 1.4319100 -0.408328595
由 reprex package (v2.0.0)
于 2021-07-08 创建
我正在管理包含 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) 模型计算的，滚动 window 10 个月，结束于 n.The 月 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.
如果我的 objective 是使用前 10 个 obs 计算 11st 的平均值，我可以简单地调用以下内容：
Mean= slide_index_dbl(Y, date, mean, .before = months(10), .after = months(-1), .complete = T)
然而，在这种情况下，objective是运行一个AR模型，使用所有之前的10个obs并使用这个估计模型来预测第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(actual data)和Predicted Y的差值。话虽如此，问题是如何得到Predicted Y。在我们预测Y之前，我们需要先训练一个AR模型，通过feed previous 10个月的数据。一旦你得到这个模型，你就可以用这个模型预测11st obs，那就是Predicted Y。最后的尝试是计算它与实际Y之间的差异，这称为shock。
数值示例如 follows:in 为了获得 2012-11-01 的输出 0.074，我需要使用从 2012-01-01 到2012-10-01，即 Y 从 0.15 到 1.85，N 从 1 到 10。训练完这个模型后，我用它来预测 2012-11-01 的新值 (0.2)。最终输出“Shock”是 2012-11-01 实际 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 与预测值 (3.46-4.1=-0.65) 之间的差异。
我不知道如何编写这样的函数(AR_2)并在slide_index_dbl中调用它。请记住 AR_2 中有一个趋势项 a3*Tt，我不确定如何对其建模。
以下是我试过的。我使用 lm 来实现 AR 模型而不是 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）。不像前面的例子，我们只是计算一个单一变量的滚动平均值，在这种情况下，为了运行这样的回归，我们需要在每个滚动window中有四个变量，但是滑块只有一个变量输入.即使我们编辑“提取”函数将四个变量减少为两个（lag_1 和 lag_2 可以写成 lag(Y,1) 和 lag(Y,2)，我们仍然有 Y 和 N两个变量输入)
关于趋势项，据我了解，例如2012-11-01，为了运行回归，我需要前10分钟N，即1,2,3。 ...10，即obs的序号。对于 2012-12-01，它还应该包括来自 1,2,3...10 的 N，但是，根据我的代码，它是 2,3,4...11。但是添加常数 1（2,3...11 对 1,2,...10）不会影响回归系数，对吧？老实说，我不太确定这种情况下的趋势术语，如果您了解原作者，请随时更改它。
你可以这样做。您的滞后 2 [A(2) 模型] 问题将由 order = c(2,0,0) 解决，您不必明确地这样做。 arima() 函数与 predict().
library(tidyverse)
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)
library(lubridate, warn.conflicts = F)
library(slider)
data %>%
mutate(pred_roll_10 = slide_index_dbl(Y, .i = date,
~ predict(arima(.x, c(2,0,0), method = 'ML'),1)$pred[1],
.before = months(10),
.after = months(-1),
.complete = T),
shock = Y - pred_roll_10)
#> Warning in log(s2): NaNs produced
#> Warning in log(s2): NaNs produced
#> Warning in log(s2): NaNs produced
#> Warning in log(s2): NaNs produced
#> date Y pred_roll_10 shock
#> 1 2012-01-01 -0.5604756 NA NA
#> 2 2012-02-01 -0.7906531 NA NA
#> 3 2012-03-01 0.7680552 NA NA
#> 4 2012-04-01 0.8385636 NA NA
#> 5 2012-05-01 0.9678513 NA NA
#> 6 2012-06-01 2.6829163 NA NA
#> 7 2012-07-01 3.1438325 NA NA
#> 8 2012-08-01 1.8787713 NA NA
#> 9 2012-09-01 1.1919184 NA NA
#> 10 2012-10-01 0.7462564 NA NA
#> 11 2012-11-01 1.9703382 0.6951079 1.275230339
#> 12 2012-12-01 2.3301521 2.0524671 0.277684986
#> 13 2013-01-01 2.7309235 1.9333295 0.797594017
#> 14 2013-02-01 2.8416062 2.0486834 0.792922863
#> 15 2013-03-01 2.2857651 1.9959595 0.289805588
#> 16 2013-04-01 4.0726782 1.7514948 2.321183420
#> 17 2013-05-01 4.5705287 3.6469980 0.923530738
#> 18 2013-06-01 2.6039116 4.0585260 -1.454614411
#> 19 2013-07-01 3.3052675 1.6480334 1.657234048
#> 20 2013-08-01 2.8324760 2.9543635 -0.121887466
#> 21 2013-09-01 1.7646523 2.8739943 -1.109341911
#> 22 2013-10-01 1.5466774 2.7845341 -1.237856702
#> 23 2013-11-01 0.5206730 2.5894221 -2.068749114
#> 24 2013-12-01 -0.2082182 1.3652172 -1.573435497
#> 25 2014-01-01 -0.8332575 0.3654956 -1.198753080
#> 26 2014-02-01 -2.5199508 -0.5625776 -1.957373251
#> 27 2014-03-01 -1.6821638 -2.6388755 0.956711703
#> 28 2014-04-01 -1.5287907 -0.7131115 -0.815679154
#> 29 2014-05-01 -2.6669276 -1.2140023 -1.452925253
#> 30 2014-06-01 -1.4131127 -2.6421098 1.228997135
#> 31 2014-07-01 -0.9866485 -1.0140132 0.027364706
#> 32 2014-08-01 -1.2817199 -0.8209904 -0.460729511
#> 33 2014-09-01 -0.3865943 -1.2655968 0.879002504
#> 34 2014-10-01 0.4915392 -1.1977841 1.689323273
#> 35 2014-11-01 1.3131203 -0.4463373 1.759457622
#> 36 2014-12-01 2.0017605 0.9931931 1.008567426
#> 37 2015-01-01 2.5556782 1.7819117 0.773766509
#> 38 2015-02-01 2.4937665 2.3699293 0.123837184
#> 39 2015-03-01 2.1878038 2.1804620 0.007341787
#> 40 2015-04-01 1.8073328 1.5452418 0.262091026
#> 41 2015-05-01 1.1126258 1.2644697 -0.151843849
#> 42 2015-06-01 0.9047086 0.2974574 0.607251157
#> 43 2015-07-01 -0.3606878 0.7761950 -1.136882824
#> 44 2015-08-01 1.8082682 -1.1719233 2.980191448
#> 45 2015-09-01 3.0162302 1.9584189 1.057811310
#> 46 2015-10-01 1.8931216 2.5356767 -0.642555104
#> 47 2015-11-01 1.4902368 1.2115509 0.278685838
#> 48 2015-12-01 1.0235814 1.4319100 -0.408328595
由 reprex package (v2.0.0)
于 2021-07-08 创建