[R] State space AR models in R: some examples
Pablo Almaraz
palmaraz at gmail.com
Thu Apr 27 16:55:19 CEST 2006
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
More information about the R-help
mailing list