ar.yw如何估計方差

在R中，函數ar.yw如何估計方差？具體來說，數字「var.pred」來自哪裏？它似乎不是來自通常的YW估計的方差，也不是平方殘差的總和除以df（即使對於df應該是什麼意見不一致，沒有一個選項給出了等同於var.pred的答案）。是的，我知道有比YW更好的方法;試圖弄清楚R在做什麼。ar.yw如何估計方差

set.seed(82346) 
temp <- arima.sim(n=10, list(ar = 0.5), sd=1) 
fit <- ar(temp, method = "yule-walker", demean = FALSE, aic=FALSE, order.max=1) 

## R's estimate of the sigma squared 
fit$var.pred 
## YW estimate 
sum(temp^2)/10 - fit$ar*sum(temp[2:10]*temp[1:9])/10 
## YW if there was a mean 
sum((temp-mean(temp))^2)/10 - fit$ar*sum((temp[2:10]-mean(temp))*(temp[1:9]-mean(temp)))/10 
## estimate based on residuals, different possible df. 
sum(na.omit(fit$resid^2))/10 
sum(na.omit(fit$resid^2))/9 
sum(na.omit(fit$resid^2))/8 
sum(na.omit(fit$resid^2))/7

來源

2016-11-18 chucky

A guess ...這是預測包中的函數嗎？你知道R包都是開源的嗎？ –

@ 42，不，這是「stats」包的一部分，即給我們「lm」和所有其他基本統計功能的相同包。因此，人們可能會認爲，即使它是開源的，它會做一些合理的事情。 – chucky

如果沒有記錄，需要閱讀代碼。

?ar.yw

其中說：「在ar.yw創新的方差矩陣是從擬合係數和x的自協方差計算出來的。如果這還不夠的解釋，那麼你需要看一下代碼：

methods(ar.yw) 
#[1] ar.yw.default* ar.yw.mts*  
#see '?methods' for accessing help and source code 

getAnywhere(ar.yw.default) 
# there are two cases that I see 
x <- as.matrix(x) 
nser <- ncol(x) 
if (nser > 1L) # .... not your situation 
#.... 
else{ 
    r <- as.double(drop(xacf)) 
    z <- .Fortran(C_eureka, as.integer(order.max), r, r, 
     coefs = double(order.max^2), vars = double(order.max), 
     double(order.max)) 
    coefs <- matrix(z$coefs, order.max, order.max) 
    partialacf <- array(diag(coefs), dim = c(order.max, 1L, 
     1L)) 
    var.pred <- c(r[1L], z$vars) 
    #....... 
    order <- if (aic) 
     (0L:order.max)[xaic == 0L] 
    else order.max 
    ar <- if (order) 
     coefs[order, seq_len(order)] 
    else numeric() 
    var.pred <- var.pred[order + 1L] 
    var.pred <- var.pred * n.used/(n.used - (order + 1L))

所以你現在需要找到Fortran代碼爲C_eureka。我想我在這裏找到它：https://svn.r-project.org/R/trunk/src/library/stats/src/eureka.f這是我認爲返回var.pred估計值的代碼。我不是一個時間系列的人，你有責任回顧這個過程，以適應你的問題。

 subroutine eureka (lr,r,g,f,var,a) 
c 
c  solves Toeplitz matrix equation toep(r)f=g(1+.) 
c  by Levinson's algorithm 
c  a is a workspace of size lr, the number 
c  of equations 
c 
    snipped 
c estimate the innovations variance 
     var(l) = var(l-1) * (1 - f(l,l)*f(l,l)) 
     if (l .eq. lr) return 
     d = 0.0d0 
     q = 0.0d0 
     do 50 i = 1, l 
      k = l-i+2 
      d = d + a(i)*r(k) 
      q = q + f(l,i)*r(k) 
    50  continue

來源

2016-11-18 17:36:56

ar.yw如何估計方差

回答

相關問題