[R] AR(2) modelling

Christophe Dutang dutangc at gmail.com
Fri Nov 13 23:30:14 CET 2009


Hi useRs,

I'm trying to fit a basic AR(2) model with the 'ar' function. And when  
I try to check the value of the coefficients, I could not find the  
same value as the 'ar' function.

Here is my example:
myserie <- c(212, 205, 210, 213, 217, 222, 216, 218, 220, 212, 215, 236)

#plot(myserie, type="l")

myserieminus0 <- tail(myserie, -2)
myserieminus1 <- tail(head(myserie, -1), -1)
myserieminus2 <- head(myserie, -2)


###Yule Walker equations

r1 <- cor(myserieminus0, myserieminus1)
r2 <- cor(myserieminus0, myserieminus2)

#method 1
phihat1 <- r1*(1-r2)/(1-r1^2)
phihat2 <- (r2-r1^2)/(1-r1^2)

#method 2
bigR <- cbind(c(1, r1), c(r1, 1))
smallr <- c(r1, r2)
ressolve <- solve(bigR, smallr)

resaryw <- ar(myserie, method="yule-walker", order.max=2, aic=FALSE)

cat("\t\tmanual YW 1\t\tmanual YW 2\t\tar YW\n")
cat("first coef", phihat1,"\t", ressolve[1],"\t\t", resaryw$ar[1], "\n")
cat("first coef", phihat2,"\t", ressolve[2],"\t", resaryw$ar[2], "\n\n")

 >		manual YW 1		manual YW 2		ar YW
 > first coef 0.2941808 	 0.2941808 		 0.1869641
 > first coef -0.1316839 	 -0.1316839 	 -0.1038001

I do not understand why the "yule-walker" does not solve exactly the  
Yule-Walker equations. A reference can be found here http://www-stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/YWSourceFiles/YW-Eshel.pdf 
  .

In the R source code (<src>/src/library/stats/R), the file ar.R  
contains the ar.yw.default function implements the function. For  
univariate case (line 130), r_eureka function is used, which seems to  
be implemented in the eureka.f function.

       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

is supposed to implement the Yule-Walker equations...

Any help is welcome.

Just to be sure, I can do something I try to reconcile the ordinary  
least square methods. And it works!

Christophe


PS : OLS code

###Ordinary Least Square

reslm <- lm(myserieminus0 ~ myserieminus1 + myserieminus2)
#summary(reslm)
coef1ols <- reslm$coefficients["myserieminus1"]
coef2ols <- reslm$coefficients["myserieminus2"]

resarols <- ar(myserie, method="ols", order.max=2, aic=FALSE)
cat("\t\tmanual ols\t\tar ols\n")
cat("first coef", coef1ols,"\t", resarols$ar[1], "\n")
cat("first coef", coef2ols,"\t", resarols$ar[2], "\n\n")






--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr




More information about the R-help mailing list