修改通常的hpfilter函数以忽略na值。

Question

我是一名新的R语言用户，想要快速学习，但是我自己无法解决这个问题。我的主要工作涉及经济时间序列 - 因此，我尝试将我的数据集维护在xts多列格式中，例如：

> head(USDATAq)
         tq   ngdp    rgdp  profit
1947 Q1   0  237.2  1770.7    20.7
1947 Q2   1  240.4  1768.0    23.9
1947 Q3   2  244.5  1766.5    23.8
1947 Q4   3  254.3  1793.3    25.5
1948 Q1   4  260.3  1821.8    29.4
1948 Q2   5  267.3  1855.3    31.2

我使用hpfilter函数进行过滤。在网站的其他地方，我找到了这个实现，它使用coredata函数将hpfilter应用于xts对象：

hpfilter <- function(x, lambda=2){
  eye <- diag(length(x))
  dcrossprod <- crossprod(diff(eye, lag=1, d=2))
  coredata(x) <- solve(eye + lambda * dcrossprod, coredata(x))
  return(x)
}

我的问题是：

如何修改函数，以使其与具有NA观测值的变量一起工作（目前，如果存在任何NA，则它会计算整个日期范围的NA）？

我可以将数据集传递为na.omit（USDATAq），这有效，但这会限制数据集中的所有变量到最小观测量。 但是，不同的变量可用直到不同的日期，然后是NA。 我最终想要在循环或mapply中将函数应用于数据集的每一列，以便该函数使用该系列的所有可用观测结果返回每个已过滤系列。

- Tatha

---

乍一看（虽然不熟悉此领域，但能够阅读R代码），似乎这既是一个统计或定义问题，也是一个R编码问题。您希望依赖滞后变量的方法如何工作？是否应该删除所有在缺失值的一个单位内的行[这个猜测基于对http://en.wikipedia.org/wiki/Hodrick%E2%80%93Prescott_filter的快速浏览]？这样做是否合理，或者会因为某些形式的缺失而导致问题...?? - Ben Bolker

猜测你是在提到我之前关于修改Farnsworth的初始hpfilter以保留ts / zoo信息的讨论，我现在刚刚给出了链接。 - Matt Bannert

是的ran2，我接受了你的修改，将coredata包含在定义中。谢谢你的提示。 - Tatha

谢谢Ben。是的，你说得对——hpfilter使用历史数据；因此，缺失的观测值在逻辑上比较棘手。但是，在这里，我只是指开始和结束时的缺失数据（也就是目前可用数据的范围而已）。中间的缺失观测值将是完全不同的问题。通常，这不是主流宏观经济数据的问题。因此，我希望该函数能够返回原始数据范围允许的日期的过滤数据。 - Tatha

Answer 1

感谢@ran2的建议。我按照您的建议进行了操作，并设法解决了问题，但是方式相当不优雅。首先，我无法让任何“apply family”函数在xts对象上正确工作，并保持其结构。使用应用(x，MARGIN=2，..)进行逐列应用的纯净apply表现出了希望，但在“coredata”语句处停滞不前。 lapply等产生了混乱的列表。

然后我转向for循环。但由于x<-na.omit(x)会改变变量的长度，因此它无法在循环内替换原始变量。

> for(i in 1:ncol(USDATAq)) {
+ USDATAq[,i]<-hpfilter(USDATAq[,i])
+ }

NextMethod(.Generic)中存在错误： 需要替换的项目数量不是多个替换长度的倍数。

因此，我不得不向hpfilter添加不雅的代码，以将结果与原始数据“合并”（包括NA），然后返回该变量。这种合并通过日期（因此是长度）来匹配2个变量，并将NA填充到结果中。然后，这个结果可以在循环中替换原始数据。总之，我必须修改hpfilter以：

hpfilter <- function(x,lambda=2){
y<-na.omit(x)
eye <- diag(length(y))
coredata(y) <- solve(eye + lambda * crossprod(diff(eye, lag=1, d=2)), coredata(y))
xy<-merge(x,y) 
return(xy[,2])
}

然后使用上述循环，最终得到无错误的结果。然而，我的R语言知识非常基础，可能会有更简单的方法。但至少我现在可以继续了。感谢所有人指导我走向正确的方向。如果对以上代码仍有改进意见，我仍然欢迎。

- Tatha

---

在编程中，保持良好的精神状态并利用提示来解决自己的问题是非常重要的。 - Ben Bolker

谢谢Ben。一切都好，暂时告一段落。还要感谢你上面提到的重新对齐技巧。我正好卡在那里了……在某个阶段，核心R将受益于为时间序列提供一个特殊的“日历”工作空间，其中将自动同步变量到日历日期（一旦定义）。 - Tatha

Answer 2

我认为你的想法是对的。为什么不在这个函数内部添加 na.omit 呢？就在创建 eye 矩阵之前加上 x<-na.omit(x)。然后，只需向其传递一元系列即可，而不是整个数据框。换句话说: 将函数保持不变，添加 na.omit 并将其与 lapply (或适合你的任何形式的 apply 家族 (sapply,tapply,lapply) 结合起来。

- Matt Bannert

---

唯一棘手的部分可能是重新填充NA，以便在之后正确匹配列...即类似于filtered.x <- x; filtered.x[!is.na(x)] <- hpfilter(x[!is.na(x)],...)这样的内容。 - Ben Bolker

真的。但这取决于你想做什么。有时我更喜欢一个单变量序列的list，而不是一个data.frame。特别是在绘图方面。然而，在动态线性模型中，例如使用dynlm，后者更方便。不过，很高兴看到OP解决了问题。 - Matt Bannert

@hans0l0 如果我们将NA替换为0会怎么样？如果需要的话，很容易再将0替换回NA。时间序列有时会用0填充以避免端点问题。 - Anusha

Answer 3

使用zoo对象时，使用attributes()比使用coredata()稍微更加清晰，然后您可以直接合并回zoo对象中。（我尚未尝试过xts对象的情况）：

hpfilter <- function(x,lambda=1600){
    y<-na.omit(x)
    eye <- diag(length(y))
    result <- solve(eye+lambda*crossprod(diff(eye,lag=1,d=2)),y)
    attributes(result) <- attributes(y)
    return(result)
}

- user2750268