Robert J. Kindman
2014-Mar-06 15:56 UTC
[R] Errors Calculating MVN Likelihood of Time Series with AR(1) Errors
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. Thank you very much, Robert -- Robert Kindman Harvard College, Class of 2014 rkindman@college.harvard.edu +1.919.599.1921 [[alternative HTML version deleted]]
peter dalgaard
2014-Mar-06 18:16 UTC
[R] Errors Calculating MVN Likelihood of Time Series with AR(1) Errors
On 06 Mar 2014, at 16:56 , Robert J. Kindman <RKindman at 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 > rkindman at college.harvard.edu > +1.919.599.1921 > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at 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.mes at cbs.dk Priv: PDalgd at gmail.com