[R] State space AR models in R: some examples

[Gabor Grothendieck]

Check out the sspir package and
http://www.jstatsoft.org/index.php?vol=16


On 4/27/06, Pablo Almaraz <palmaraz at gmail.com> wrote:
> Hi all,
>
> Does anyone have an example of an autoregressive (AR) time-series model
> specified as a state space model in R? That is, I want to go beyond the
> locally linear (constant) model, and fit the following Gaussian AR state
> process model:
>
> Xt = a + (1+b)*Xt-1 + epsilon
>
> ,where the model for the observation process is
>
> Yt = Xt + tau
>
> I have information of the tau's (observation variance) for each
> observation in the time-series, and it would be perfect to include this
> information during the fitting routine. I have actually coded this as a
> WinBUGS code (pasted below), but I'm not quite sure it works as it
> should. I would be extremely thanked if anyone could submit an example
> of an R code fitting the above problem. Gibb's sampler for solving the
> problem would be great, I'm not sure whether the Kalman filter would
> work well with only 30 data points (?). Additional details, corrections
> and/or help would probably save my life at least for a while.
>
> Thank you all
> Cheers
>
> Pablo
>
> WinBUGS state-space R code:
> ######
> model;
>    {
>
> # Parameters and priors
>        alpha ~ dnorm(0,0.000001)   # Intrinsic rate of increase
>        b ~ dnorm(0,0.000001)
>        beta[1] <- b-1   # First-order density-dependence
>        sigma ~ dunif(0, 1000)   # State process SD
>        isigma2 <- pow(sigma, -2)   # State process 1/var
>        # isigma2 ~ dgamma(0.01,0.000001)
>
> # Initial state value
>        n.exp[1] ~ dnorm(n[1],tau[1])
>
> # State process model
>        for(j in 1:(N-1)){
>            n.exp.mu[j+1] <- alpha + b*n.exp[j]  # First-order Gompertz
> model
>            n.exp[j+1] ~ dnorm(n.exp.mu[j+1], isigma2)
>        }
>
> # Observation process model
>        for(j in 1:(N-1)){
>            n[j+1] ~ dnorm(n.exp[j+1],tau[j+1])
>        }
>    }
>
> # Loge-transformed and standardized time-series data
>
> list(N=28,
> n=c(-0.24645, 0.015312, 0.442262, -0.05879, -0.17308, -0.03778,
> 0.120961, -0.04383, 0.002507, 0.073278, -0.11684, 0.003657, -0.07375,
> 0.05006, -0.04489, -0.00826, -0.06713, 0.682228, 0.032058, -0.33254,
> -0.50432, 0.176914, 0.249793, 0.01672, -0.30581, -0.19617, 0.158579,
> 0.185296),
> tau=c(2.38351, 2.351379, 49.12811, 10.01703, 11.68982, 3.846619,
> 1.999254, 1.6685, 3.011932, 5.661051, 168.2524, 1.581, 25.74985,
> 50.29332, 3.03117, 7.65013, 3.376606, 17.34871, 4.215985, 2.455294,
> 7.685724, 1.918054, 5.588953, 8.503541, 0.5666, 0.923611, 4.986243,
> 10.36613))
>
> # Inits for MC 1
> list(alpha = 0.5, b = -1, sigma = 0.5)
>
> # Inits for MC 2
> list(alpha = 1, b = 0.01, sigma = 1)
>
> # Inits for MC 3
> list(alpha = 0.01, b = 1, sigma = 0.01)
>
> ### End (not run)
>
> --
>
> Pablo Almaraz García
> Estación Biológica de Doñana (CSIC)
> Pabellón del Perú, Avda. Mª Luísa s/n
> E-41013, Sevilla
> SPAIN
>
> E-mail: almaraz[AT]ebd[DOT]csic[DOT]es
> webpage: http://www.almaraz.org
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>