R2WinBUGS: Autoregressive model of order 1, AR(1)

Let us assume that y₁,…,y_N are observed and that a first order autoregressive model, AR(1), is entertained

y_t ~ N(a+by_t-1,t^-2)

for t=1,…,N. Additionally, two different prior specifications will be examined:

π(a,b,t²,y₀) = π(a)π(b)π(t²)π(y₀)

PRIOR 1: Conditionally conjugate components

a ~ N(0,1000000)

b ~ N(0,1000000)

t² ~ G(0.01,0.01)

y₀ ~ N(0,1000000)

PRIOR 2: Stationarity restriction

a ~ N(0,1000000)

b ~ N(0,1000000)1(0<b<1)

t² ~ G(0.01,0.01)

y₀ ~ N(a/(1-b), t^-2/(1-b²))

WINBUGS CODE FOR THE MODEL AND PRIOR 1

# This model is stored in file “AR1-1.bug”

model{

# AR(1) model

mu[1] <- alpha + beta*y0

y[1] ~ dnorm(mu[1],tau2)

for(t in 2:N) {

mu[t] <- alpha + beta*y[t-1]

y[t] ~ dnorm(mu[t],tau2)

}

# Prior distribution

alpha ~ dnorm(0.0, 1.0E-6)

beta ~ dnorm(0.0, 1.0E-6)

tau2 ~ dgamma(0.01, 0.01)

y0 ~ dnorm(0.0,1.0E-6)

}

WINBUGS CODE FOR THE MODEL AND PRIOR 2

# This model is stored in file “AR1-2.bug”

model{

# AR(1) model

mu[1] <- alpha + beta*y0

y[1] ~ dnorm(mu[1],tau2)

for(t in 2:N) {

mu[t] <- alpha + beta*y[t-1]

y[t] ~ dnorm(mu[t],tau2)

}

# Prior distribution

m <- alpha/(1-beta)

p <- (1-beta*beta)*tau2

alpha ~ dnorm(0.0, 1.0E-6)

beta ~ dnorm(0.0, 1.0E-6)I(0.0,1.0)

tau2 ~ dgamma(0.01, 0.01)

y0 ~ dnorm(m,p)

}

Evoking WinBUGS from R

# Loading the R2WingBUGS library

library(R2WinBUGS)

# Simulating N=100 observations form an AR(1) process

set.seed(1245)

N = 100

a = 1.00

b = 0.95

s2 = 1.00

y0 = rnorm(1,a/(1-b),sqrt(s2/(1-b^2)))

y = rep(0,N)

y[1] = rnorm(1,a+b*y0,sqrt(s2))

for (t in 2:N)

y[t] = rnorm(1,a+b*y[t-1],sqrt(s2))

ts.plot(y)

# Posterior inference: the conditionally conjugate case

# -----------------------------------------------------

data <- list("y","N")

inits <- function(){

list(alpha=1.0,beta=0.95,tau2=1.0,y0=20)

}

# Running WinBUGS

AR11 = bugs(data,inits,model.file="AR1-1.bug",

parameters=c("alpha","beta","tau2","y0"),

n.chains=1,n.iter=10000,n.burnin=5000,n.thin=1,

bugs.directory="c:/Program Files/WinBUGS14/",codaPkg=FALSE)

# See what the object AR11 contains

AR11

# See the list of components in AR11

names(AR11)

# Posterior summary

names = c("alpha","beta","tau2")

true = c(a,b,1/s2)

par(mfrow=c(3,3))

for (i in 1:3){

ts.plot(AR11$sims.array[,1,i],xlab="iterations",ylab="",main=names[i])

abline(h=true[i],col=2,lwd=3)

acf(AR11$sims.array[,1,i],xlab="lag",ylab="ACF",main="")

hist(AR11$sims.array[,1,i],xlab="",ylab="",main="",prob=T)

abline(v=true[i],col=2,lwd=3)

}

# Posterior inference: the stationary case

# ----------------------------------------

data <- list("y","N")

inits <- function(){

list(alpha=1.0,beta=0.95,tau2=1.0,y0=20)

}

AR12 = bugs(data,inits,model.file="AR1-2.bug",

parameters=c("alpha","beta","tau2","y0"),

n.chains=1,n.iter=255000,n.burnin=5000,n.thin=50,

bugs.directory="c:/Program Files/WinBUGS14/",codaPkg=FALSE)

# Posterior summary

names = c("alpha","beta","tau2")

true = c(a,b,1/s2)

par(mfrow=c(3,3))

for (i in 1:3){

ts.plot(AR12$sims.array[,1,i],xlab="iterations",ylab="",main=names[i])

abline(h=true[i],col=2,lwd=3)

acf(AR12$sims.array[,1,i],xlab="lag",ylab="ACF",main="")

hist(AR12 $sims.array[,1,i],xlab="",ylab="",main="",prob=T)

abline(v=true[i],col=2,lwd=3)

}

# Comparing posterior distributions based on both prior distributions

par(mfrow=c(2,2))

boxplot(AR11$sims.array[,1,1],AR12$sims.array[,1,1],

names=c("Conjugate","Stationary"),main="alpha")

boxplot(AR11$sims.array[,1,2],AR12$sims.array[,1,2],

names=c("Conjugate","Stationary"),main="beta")

boxplot(AR11$sims.array[,1,3],AR12$sims.array[,1,3],

names=c("Conjugate","Stationary"),main="tau2")

boxplot(AR11$sims.array[,1,4],AR12$sims.array[,1,4],

names=c("Conjugate","Stationary"),main="y0")