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)