bad estimate of error variance for AR(1)

[Timothée Flutre]

Hello,
I am playing with the AR(1), but don't quite understand how to use inla for it:

• for all t in {2,...,T}, x_t = mu + rho x_{t-1} + e_t with e_t ~ N(0, sigma^2)
  
• x_1 ~ N(mu, sigma^2 / (1 - rho^2))

I am simulating data according to this model (R version 3.3.1):

set.seed(1)

T <- 10^3

mu <- 20

rho <- 0.7

sigma2 <- 5

x <- mu + arima.sim(model=list(ar=rho), n=T, sd=sqrt(sigma2))

And I am fitting the same model to these data, both with arima and inla (INLA version 0.0-1468872408, dated 2016-07-18):

fit1 <- arima(x=x, order=c(1,0,0))

fit2 <- inla(formula=x ~ 1 + f(i, model="ar1"), data=data.frame(x=x, i=1:length(x)))

Given the fairly large sample size, I expect the estimates of mu, rho and sigma to be quite accurate.

Indeed, the results are fine for arima:

Coefficients:

      ar1     intercept

      0.6511  19.8572

s.e.  0.0240   0.2093

sigma^2 estimated as 5.353: log likelihood = -2258.03, aic = 4522.07

Here are the results for inla:

            mean    sd     0.025quant 0.5quant 0.975quant mode    kld

(Intercept) 19.8572 0.2106 19.4426    19.8572  20.2711    19.8572 0

                                        mean       sd         0.025quant  0.5quant   0.975quant  mode

Precision for the Gaussian observations 1.982e+04  2.033e+04  1412.0635   1.378e+04  7.340e+04   3876.447

Precision for i                         1.072e-01  7.500e-03  0.0929      1.071e-01  1.224e-01       0.107

Rho for i                               6.528e-01  2.390e-02  0.6052      6.530e-01  6.989e-01       0.653

The estimates of mu and rho are fine.

1. However, the "precision for i" does correspond to 1 / sigma^2, right?
2. If yes, why isn't the true value included in the 95% credible interval [1/0.1224=8.17 ; 1/0.0929=10.76]?
   
3. And what does the "precision for the Gaussian observations" correspond to?

My guess is that the way I am currently using inla assumes a latent AR(1) process (as written in the beginning), to which is superimposed some random errors, as if we were observing y_t = x_t + w_t where w_t ~ N(0, sigma_w^2) and the "precision for the Gaussian observations" corresponds to 1/sigma_w^2. As I simulate data with sigma_w^2 = 0, it's fine if "precision of the Gaussian observations" is very large.

But then, why is the posterior of sigma^2 so off-target? Does it mean I should tune the hyper-priors of theta1 = log(kappa), where kappa is the marginal precision, to flatten the prior on kappa?

And is there a way with inla to directly fit the AR(1) as the observed process instead of as a latent one, e.g. by fixing sigma_w^2 to 0?

Thanks in advance,

Tim

[Elias T. Krainski]

Two points:

1. you have to fix the precision for the likelihood on a big value 

 inla(..., control.family=list(list(hyper=list(theta=list(initial=10, fixed=TRUE)))))

2. ar1 in INLA is parametrized as marginal variance, see inla.doc('ar1')

Elias

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.

[Håvard Rue]

On Thu, 2016-10-13 at 22:49 -0300, Elias T. Krainski wrote:
> Two points:
> 1. you have to fix the precision for the likelihood on a big value 
>  inla(..., control.family=list(list(hyper=list(theta=list(initial=10,
> fixed=TRUE)))))
> 2. ar1 in INLA is parametrized as marginal variance, see
> inla.doc('ar1')

yes, here is an example

> n = 100
> idx = 1:n
> y = arima.sim(n, model = list(ar = 0.9))
> 1/var(y)
[1] 0.4154346454
> r = inla(y ~ 1 + f(idx, model="ar1"), data = data.frame(y,idx),
control.family = list(hyper = list(prec = list(initial = 10,
fixed=TRUE))))
> summary(r)
::::
Model hyperparameters:
                    mean     sd 
Precision for idx 0.3968 0.1138 
Rho for idx       0.8188 0.0519 

--
Håvard Rue
Department of Mathematical Sciences
Norwegian University of Science and Technology
N-7491 Trondheim, Norway
Voice: +47-7359-3533 URL : http://www.math.ntnu.no/~hrue
Mobile: +47-9260-0021 Email: havar...@math.ntnu.no

R-INLA: www.r-inla.org

[Timothée Flutre]

Thanks! I think I get it now.