The data are the HadCRUT3v data set, which I actually took from this blogpost. Here I reworked the analysis for the Bayesian setting. The data decribe mean global annual temperature anomalies from 1850 to 2013, provided at yearly intervals.

# 1 The data

   gtemp <- read.csv("http://www.petrkeil.com/wp-content/uploads/2016/01/gtemp.csv")
plot(Annual~Year, data=gtemp, type="b")

Letâ€™s explore the temporal autocorrelation a bit:

   par(mfrow=c(1,2))
plot(gtemp$Annual[-164], gtemp$Annual[-1],
xlab="T in Year", ylab="T in Year-1")
acf(gtemp$Annual, main="") # 2 First-order autoregressive (AR1) model in JAGS Formal definition of the 1st-order autoregressive (AR1) model: $y_t = \beta_0 + \beta_1 y_{t-1} + \epsilon_i$ ## 2.1 Prepare the data  ar.data <- list(Annual = gtemp$Annual,
N = nrow(gtemp))

## 2.2 The model

We will use the R2jags library:

   library(R2jags)

JAGS model specification:

cat("
model
{
# priors
tau <- 1/(sigma*sigma)
sigma ~ dunif(0, 100)
beta0 ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
mu[1] ~ dnorm(0, 0.001)

# likelihood
eps[1] <- Annual[1] - mu[1]

for(t in 2:N)
{
mu[t] <- beta0 + beta1*Annual[t-1]
eps[t] <- Annual[t] - mu[t]
}
for(t in 1:N)
{
Annual[t] ~ dnorm(mu[t], tau)
}
}
", file="ar1.txt")

Fitting the model:

   model.fit <- jags(data=ar.data,
model.file="ar1.txt",
parameters.to.save=c("beta0", "beta1"),
n.chains=3,
n.iter=2000,
n.burnin=1000)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 164
##    Unobserved stochastic nodes: 4
##    Total graph size: 644
##
## Initializing model
   model.fit
## Inference for Bugs model at "ar1.txt", fit using jags,
##  3 chains, each with 2000 iterations (first 1000 discarded)
##  n.sims = 3000 iterations saved
##           mu.vect sd.vect     2.5%      25%      50%      75%    97.5%
## beta0      -0.004   0.009   -0.022   -0.010   -0.004    0.002    0.014
## beta1       0.940   0.030    0.880    0.920    0.941    0.961    1.000
## deviance -288.841   2.908 -292.465 -290.991 -289.511 -287.356 -281.623
##           Rhat n.eff
## beta0    1.001  3000
## beta1    1.001  3000
## deviance 1.001  3000
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 4.2 and DIC = -284.6
## DIC is an estimate of expected predictive error (lower deviance is better).
   mus <- jags(data=ar.data,
model.file="ar1.txt",
parameters.to.save=c("mu"),
n.chains=3,
n.iter=2000,
n.burnin=1000,
DIC=FALSE)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 164
##    Unobserved stochastic nodes: 4
##    Total graph size: 644
##
## Initializing model
   mus <- data.frame(mus$BUGSoutput$summary)

plot(Annual~Year, data=gtemp, type="b")
lines(gtemp$Year, mus$X50., col="red", lwd=2)