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

 

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

 

yt  ~ N(a+byt-1,t-2)

 

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

 

                                                π(a,b,t2,y0)  = π(a)π(b)π(t2)π(y0)

 

PRIOR 1: Conditionally conjugate components

 

                                                a   ~ N(0,1000000)

                                                b   ~ N(0,1000000)

                                                t2  ~ G(0.01,0.01)

                                                y0   ~ N(0,1000000)

 

PRIOR 2: Stationarity restriction

 

                                                a   ~ N(0,1000000)

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

                                                t2    ~ G(0.01,0.01)

                                                y0  ~ N(a/(1-b), t-2/(1-b2))

 

 

WINBUGS CODE FOR THE MODEL AND PRIOR 1

 

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)

}

 

# Simulated data: N=100, alpha=1.00, beta=0.95, tau2=1.00 and y0= 21.34232

list(y=c(1.98194, 1.83725, 1.74383, 1.70332, 1.64980, 1.69990, 1.80716, 1.81172, 1.85085, 1.80320, 1.73829, 1.83654, 1.84023, 1.69967, 1.68530, 1.51206, 1.46801, 1.56422, 1.43848, 1.40941, 1.33075, 1.23686, 1.32284, 1.23779, 1.25595, 1.12693, 1.31803, 1.44500, 1.52974, 1.61930, 1.75544, 1.75718, 1.83342, 1.75297, 1.78100, 1.80371, 1.95254, 2.00005, 1.89442, 1.86405, 1.90428, 1.88808, 1.82147, 1.75183, 1.76600, 1.79157, 1.74878, 1.91140, 1.84219, 1.84891, 1.95664, 1.83369, 2.05875, 1.88695, 1.90774, 1.79298, 1.95589, 1.87760, 1.96444, 1.79710, 2.02845, 2.02460, 2.00059, 1.96465, 1.89335, 1.94758, 1.85429, 1.75827, 1.68618, 1.75268, 1.76941, 1.67950, 1.56109, 1.55348, 1.59788, 1.74139, 1.75932, 1.86723, 1.69596, 1.68172, 1.55348, 1.68516, 1.77895, 1.77244, 1.70491, 1.55606, 1.62143, 1.80131, 1.77279, 1.81020, 1.81793, 1.95572, 1.94176, 1.87172, 1.85153, 1.92567, 2.06492, 2.12328, 2.16769, 2.12076), N=100)

 

# Initial values

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

 

 

WINBUGS CODE FOR THE MODEL AND PRIOR 2

 

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)

}

 

# Simulated data: N=100, alpha=1.00, beta=0.95, tau2=1.00 and y0= 21.34232

list(y=c(1.98194, 1.83725, 1.74383, 1.70332, 1.64980, 1.69990, 1.80716, 1.81172, 1.85085, 1.80320, 1.73829, 1.83654, 1.84023, 1.69967, 1.68530, 1.51206, 1.46801, 1.56422, 1.43848, 1.40941, 1.33075, 1.23686, 1.32284, 1.23779, 1.25595, 1.12693, 1.31803, 1.44500, 1.52974, 1.61930, 1.75544, 1.75718, 1.83342, 1.75297, 1.78100, 1.80371, 1.95254, 2.00005, 1.89442, 1.86405, 1.90428, 1.88808, 1.82147, 1.75183, 1.76600, 1.79157, 1.74878, 1.91140, 1.84219, 1.84891, 1.95664, 1.83369, 2.05875, 1.88695, 1.90774, 1.79298, 1.95589, 1.87760, 1.96444, 1.79710, 2.02845, 2.02460, 2.00059, 1.96465, 1.89335, 1.94758, 1.85429, 1.75827, 1.68618, 1.75268, 1.76941, 1.67950, 1.56109, 1.55348, 1.59788, 1.74139, 1.75932, 1.86723, 1.69596, 1.68172, 1.55348, 1.68516, 1.77895, 1.77244, 1.70491, 1.55606, 1.62143, 1.80131, 1.77279, 1.81020, 1.81793, 1.95572, 1.94176, 1.87172, 1.85153, 1.92567, 2.06492, 2.12328, 2.16769, 2.12076), N=100)

 

# Initial values

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