Áp dụng hồi quy cửa sổ cuộn cho chuỗi XTS trong R

Question

Tôi có xts 1033 điểm trả về hàng ngày cho 5 cặp tiền tệ mà tôi muốn chạy hồi quy cửa sổ cuộn, nhưng rollapply không hoạt động cho hàm được xác định của tôi sử dụng lm(). Đây là dữ liệu của tôi:

> head(fxr) 
       USDZAR  USDEUR  USDGBP  USDCHF  USDCAD 
2007-10-18 -0.005028709 -0.0064079963 -0.003878743 -0.0099537170 -0.0006153215 
2007-10-19 -0.001544470 0.0014275520 -0.001842564 0.0023058211 -0.0111410271 
2007-10-22 0.010878027 0.0086642116 0.010599365 0.0051899551 0.0173792230 
2007-10-23 -0.022783987 -0.0075236355 -0.010804304 -0.0041668499 -0.0144788687 
2007-10-24 -0.006561223 0.0008545792 0.001024275 -0.0004261666 0.0049525483 
2007-10-25 -0.014788901 -0.0048523001 -0.001434280 -0.0050425302 -0.0046422944 

> tail(fxr) 
       USDZAR  USDEUR  USDGBP  USDCHF  USDCAD 
2012-02-10 0.018619309 0.007548205 0.005526184 0.006348533 0.0067151342 
2012-02-13 -0.006449463 -0.001055966 -0.002206810 -0.001638002 -0.0016995755 
2012-02-14 0.006320364 0.006843933 0.006605875 0.005992935 0.0007001751 
2012-02-15 -0.001666872 0.004319096 -0.001568874 0.003686840 -0.0015009759 
2012-02-16 0.006419616 -0.003401364 -0.005194817 -0.002709588 -0.0019044761 
2012-02-17 -0.004339687 -0.003675992 -0.003319899 -0.003043481 0.0000000000 

tôi có thể dễ dàng chạy một lm vào nó cho các bộ tập dữ liệu để mô hình USDZAR chống lại các cặp khác:

> lm(USDZAR ~ ., data = fxr)$coefficients 
    (Intercept)  USDEUR  USDGBP  USDCHF  USDCAD 
-1.309268e-05 5.575627e-01 1.664283e-01 -1.657206e-01 6.350490e-01 

Tuy nhiên tôi muốn chạy một cửa sổ 62 ngày lăn để có được sự phát triển của các hệ số này theo thời gian, vì vậy tôi tạo ra một dolm chức năng mà thực hiện điều này:

> dolm 
function(x) { 
    return(lm(USDZAR ~ ., data = x)$coefficients) 
} 

Tuy nhiên khi tôi chạy rollapply về vấn đề này tôi nhận được như sau:

> rollapply(fxr, 62, FUN = dolm) 
Error in terms.formula(formula, data = data) : 
    '.' in formula and no 'data' argument 

đó là mặc dù dolm (fxr) trên tác phẩm riêng của mình tốt:

> dolm(fxr) 
    (Intercept)  USDEUR  USDGBP  USDCHF  USDCAD 
-1.309268e-05 5.575627e-01 1.664283e-01 -1.657206e-01 6.350490e-01 

gì đang xảy ra ở đây? Dường như hoạt động tốt nếu dolm là hàm đơn giản hơn, ví dụ:

> dolm <- edit(dolm) 
> dolm 
function(x) { 
    return(mean(x)) 
} 
> rollapply(fxr, 62, FUN = dolm) 
        USDZAR  USDEUR  USDGBP  USDCHF  USDCAD 
2007-11-29 -1.766901e-04 -6.899297e-04 6.252596e-04 -1.155952e-03 7.021468e-04 
2007-11-30 -1.266130e-04 -6.512204e-04 7.067767e-04 -1.098413e-03 7.247315e-04 
2007-12-03 8.949942e-05 -6.406932e-04 6.637066e-04 -1.154806e-03 8.727564e-04 
2007-12-04 2.042046e-04 -5.758493e-04 5.497422e-04 -1.116308e-03 7.124593e-04 
2007-12-05 7.343586e-04 -4.899982e-04 6.161819e-04 -1.057904e-03 9.915495e-04 

Bất kỳ trợ giúp nào được đánh giá cao. Về cơ bản những gì tôi muốn là để có được trọng số cho sự hồi quy của USDZAR ~ USDEUR + USDGBP + USDCHF + USDCAD trên một cửa sổ 62 ngày lăn.

Answer 1

Có một số vấn đề ở đây:

• rollapply qua một ma trận nhưng lm đòi hỏi một data.frame.
• rollapply áp dụng chức năng cho từng cột riêng biệt trừ khi chúng tôi chỉ định by.column=FALSE.
• bạn có thể hoặc có thể không muốn kết quả là đúng phù hợp với những ngày nhưng nếu bạn sử dụng rollapplyr:

1) Kết hợp trên ta có:

dolm <- function(x) coef(lm(USDZAR ~ ., data = as.data.frame(x)))) 
rollapplyr(fxr, 62, dolm, by.column = FALSE) 

2) Cách thay thế cho số lm trong số dolm ở trên là sử dụng lm.fit làm việc trực tiếp với ma trận và cũng nhanh hơn:

dolm <- function(x) coef(lm.fit(cbind(Intercept = 1, x[,-1]), x[,1])) 

Answer 2

Tùy chọn thứ ba là cập nhật ma trận R trong phân tích QR như được thực hiện trong mã bên dưới. Bạn có thể tăng tốc độ này lên bằng cách thực hiện nó trong C++ nhưng hơn bạn sẽ cần dchud và dchdd chương trình con từ LINPACK (hoặc một chức năng để cập nhật R)

library(SamplerCompare) # for LINPACK `chdd` and `chud` 
roll_coef <- function(X, y, width){ 
    n <- nrow(X) 
    p <- ncol(X) 
    out <- matrix(NA_real_, n, p) 

    is_first <- TRUE 
    i <- width 
    while(i <= n){ 
    if(is_first){ 
     is_first <- FALSE 
     qr. <- qr(X[1:width, ]) 
     R <- qr.R(qr.) 

     # Use X^T for the rest 
     X <- t(X) 

     XtY <- drop(tcrossprod(y[1:width], X[, 1:width])) 
    } else { 
     x_new <- X[, i] 
     x_old <- X[, i - width] 

     # update R 
     R <- .Fortran(
     "dchud", R, p, p, x_new, 0., 0L, 0L, 
     0., 0., numeric(p), numeric(p), 
     PACKAGE = "SamplerCompare")[[1]] 

     # downdate R 
     R <- .Fortran(
     "dchdd", R, p, p, x_old, 0., 0L, 0L, 
     0., 0., numeric(p), numeric(p), integer(1), 
     PACKAGE = "SamplerCompare")[[1]] 

     # update XtY 
     XtY <- XtY + y[i] * x_new - y[i - width] * x_old 
    } 

    coef. <- .Internal(backsolve(R, XtY, p, TRUE, TRUE)) 
    out[i, ] <- .Internal(backsolve(R, coef., p, TRUE, FALSE)) 

    i <- i + 1 
    } 

    out 
} 

# simulate data 
set.seed(101) 
n <- 1000 
wdth = 100 
X <- matrix(rnorm(10 * n), n, 10) 
y <- drop(X %*% runif(10)) + rnorm(n) 
Z <- cbind(y, X) 

# assign other function 
dolm <- function(x) 
    coef(lm.fit(x[, -1], x[, 1])) 

# show that they yield the same 
library(zoo) 
all.equal(
    rollapply(Z, wdth, FUN = dolm, 
      by.column = FALSE, align = "right", fill = NA_real_), 
    roll_coef(X, y, wdth), 
    check.attributes = FALSE) 
#R> [1] TRUE 

# benchmark 
library(compiler) 
roll_coef <- cmpfun(roll_coef) 
dolm <- cmpfun(dolm) 
microbenchmark::microbenchmark(
    new = roll_coef(X, y, wdth), 
    prev = rollapply(Z, wdth, FUN = dolm, 
        by.column = FALSE, align = "right", fill = NA_real_), 
    times = 10) 
#R> Unit: milliseconds 
#R> expr  min   lq  mean  median   uq  max neval cld 
#R> new 8.631319 9.010579 9.808525 9.659665 9.973741 11.87083 10 a 
#R> prev 118.257128 121.734860 124.489826 122.882318 127.195410 135.21280 10 b 

Các giải pháp trên đòi hỏi bạn phải hình thành model.matrix và model.response đầu tiên nhưng đây chỉ là ba cuộc gọi (thêm một cuộc gọi đến model.frame) trước cuộc gọi đến roll_coef.