Wie schätzt ar.yw die Varianz

Question

In R, wie schätzt die Funktion ar.yw die Varianz? Wo kommt die Nummer "var.pred" konkret her? Es scheint nicht von der üblichen YW-Schätzung der Varianz zu kommen, noch von der Summe der quadrierten Residuen geteilt durch df (obwohl es Uneinigkeit darüber gibt, was das df sein sollte, gibt keine der Auswahl eine Antwort, die äquivalent zu var.pred ist). . Und ja, ich weiß, dass es bessere Methoden als YW gibt; Ich versuche nur herauszufinden, was R macht.

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 

Answer 1

Müssen Sie den Code lesen, wenn es nicht dokumentiert ist.

?ar.yw 

die sagt: „In ar.yw die Varianzmatrix der Neuerungen von den angepaßten Koeffizienten und der Autokovarianz von x berechnet wird.“ Wenn das nicht genug Erklärung ist, dann müssen Sie auf den Code suchen:

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)) 

So, jetzt Sie die Fortran-Code für C_eureka finden müssen. Ich denke, ich finde es hier: https://svn.r-project.org/R/trunk/src/library/stats/src/eureka.f Dies ist der Code, den ich denke, gibt die var.pred Schätzung zurück. Ich bin kein Zeitserien-Typ und es liegt in Ihrer Verantwortung, diesen Prozess auf Anwendbarkeit für Ihr Problem zu überprüfen.

 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