On 06 Mar 2014, at 16:56 , Robert J. Kindman <rkind...@college.harvard.edu> wrote:
> Hi all, > > I'm having trouble calculating the likelihood of a time series with AR(1) > errors. I am generating my covariance matrix according to page 2 of ( > http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-timeseries-regression.pdf), > using the library mvtnorm and the multivariate normal density function > dmvnorm(). > Here's some example code: > > library(mvtnorm) > > # Generate a basic time series with AR(1) Errors: > > t <- 1:100 > > error <- as.numeric(arima.sim(n = length(t), list(ar = c(0.8897)), sd = 10)) > > series <- 5*t + error > > # Fit the series using a basic linear model assuming errors are IID Normal > > naive.model <- lm(series ~ t -1) > > # Examine and model the residuals > > residuals <- series - t*coef(naive.model) > > residual.model <- arima(residuals, c(1,0,0), include.mean=F) > > # Construct the covariance matrix, assuming the process variance (10^2) is > known > > sigma <- diag(length(t)) > > sigma[(abs(row(sigma)-col(sigma)) == 1)] = as.numeric(coef(residual.model)) > > sigma <- sigma*10^2 > > # Calculate the MVN density... > > dmvnorm(series, t*coef(naive.model) ,sigma, log=T) > > > Without fail, I get the following error message: > > Warning message: In log(eigen(sigma, symmetric = TRUE, only.values = > TRUE)$values) : NaNs produced. > It's worth noting that the matrix from the following ( > https://stat.ethz.ch/pipermail/r-help/2007-May/131728.html) "works", but I > think is actually for an MA(1) process rather than an AR(1) process. > > I gather the message means the proposed covariance matrix may not be > invertible. This said I'm stuck on how to proceed and would be extremely > appreciative of any thoughts. > > You need to review your theory. MA-processes have band-diagonal covariance matrices, AR processes do not. The inverse of an AR covariance matrix is band diagonal, since conditional correlations are zero, but there are edge effects so that the diagonal of the inverse is not constant. Adding to that, and MA process with a neighbor correlation on the order of .9 is not possible. To paraphrase, you're doing essentially this: > M <- diag(5) > diag(M[,-1]) <- diag(M[-1,]) <- .9 > M [,1] [,2] [,3] [,4] [,5] [1,] 1.0 0.9 0.0 0.0 0.0 [2,] 0.9 1.0 0.9 0.0 0.0 [3,] 0.0 0.9 1.0 0.9 0.0 [4,] 0.0 0.0 0.9 1.0 0.9 [5,] 0.0 0.0 0.0 0.9 1.0 > eigen(M) $values [1] 2.5588457 1.9000000 1.0000000 0.1000000 -0.5588457 .... (Notice that one eigenvalue is negative, so M is not a covariance matrix. This happens already in the 3x3 case) Presumably, what you wanted is this: > M <- diag(5) > M <- .9 ^ abs(row(M)-col(M)) > M [,1] [,2] [,3] [,4] [,5] [1,] 1.0000 0.900 0.81 0.729 0.6561 [2,] 0.9000 1.000 0.90 0.810 0.7290 [3,] 0.8100 0.900 1.00 0.900 0.8100 [4,] 0.7290 0.810 0.90 1.000 0.9000 [5,] 0.6561 0.729 0.81 0.900 1.0000 > eigen(M) $values [1] 4.26213409 0.45446614 0.14592141 0.07943386 0.05804450 ..... > zapsmall(solve(M)) [,1] [,2] [,3] [,4] [,5] [1,] 5.263158 -4.736842 0.000000 0.000000 0.000000 [2,] -4.736842 9.526316 -4.736842 0.000000 0.000000 [3,] 0.000000 -4.736842 9.526316 -4.736842 0.000000 [4,] 0.000000 0.000000 -4.736842 9.526316 -4.736842 [5,] 0.000000 0.000000 0.000000 -4.736842 5.263158 > My time series knowledge hasn't really been exercised for 30 years, so I don't recall whether there is a nice formula for the entries in solve(M)... > Thank you very much, > Robert > > > -- > Robert Kindman > Harvard College, Class of 2014 > rkind...@college.harvard.edu > +1.919.599.1921 > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd....@cbs.dk Priv: pda...@gmail.com ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
