AR(1), LLM and AR(1) plus noise models

1 Global mean temperature anomaly

Global mean temperature anomaly (difference between current year and the 1951-1980 average) - Section 3.1.4 of Reich and Ghosh (2026) Bayesian Statistical Methods: With Applications to Machine Learning (2nd edition).

y  = scan("anomaly.txt")
n  = length(y)
par(mfrow=c(1,1))
plot(y,type="b",ylab="Temperature anomaly (C)",xlab="Year")

2 MODEL 1: AR(1)

Here we model data\(=\{y_1,\ldots,y_n\}\) as a Gaussian AR(1) structure. \[ y_t|y_{t-1},\rho,\sigma^2 \sim N(\rho y_{t-1},\sigma^2(1-\rho^2)), \qquad t=2,\ldots,n, \] for \(y_1 \sim N(0,\sigma^2)\). The prior for \((\rho,\sigma^2)\) is \[ (\rho,\sigma^2) \sim Beta(a,b)IG(c,d), \] for \((a,b,c,d)=(104,5,6,10)\).

prior = function(rho,sig){
  dbeta(rho,a,b)*dgamma(sig,c,d)
}

like = function(rho,sig){
  sd = sig*sqrt((1-rho^2))
  like = dnorm(y[1],0,sig)*prod(dnorm(y[2:n],rho*y[1:(n-1)],sd))
  return(like)
}

post = function(rho,sig){
  like(rho,sig)*prior(rho,sig)
}

a = 104
b = 5
c = 6
d = 10

N      = 200
sigs   = seq(0.1,1.75,length=N)
rhos   = seq(0.9,1.0,length=N)
likes  = matrix(0,N,N)
priors = matrix(0,N,N)
for (i in 1:N){
  for (j in 1:N){
    priors[i,j] = prior(rhos[i],sigs[j])
    likes[i,j] = like(rhos[i],sigs[j])
  }
}
posts = priors*likes

par(mfrow=c(1,2),mar = c(4, 4, 1, 1))
contour(rhos,sigs,priors,xlab=expression(rho),ylab=expression(sigma),drawlabels=FALSE)
contour(rhos,sigs,likes,xlab=expression(rho),ylab=expression(sigma),drawlabels=FALSE,add=TRUE,col=2)
legend("topleft",legend=c("Likelihood","Prior"),col=2:1,bty="n",lty=1)

contour(rhos,sigs,priors,xlab=expression(rho),ylab=expression(sigma),drawlabels=FALSE)
contour(rhos,sigs,posts,xlab=expression(rho),ylab=expression(sigma),drawlabels=FALSE,add=TRUE,col=2)
legend("topleft",legend=c("Posterior","Prior"),col=2:1,bty="n",lty=1)

2.1 Random-walk Metropolis-Hastings algorithm

rho.d  = 0.9
sig.d  = 1.4
sd.rho = 0.02
sd.sig = 0.05
M0     = 10000
M      = 5000
L      = 50
niter = M0+M*L
draws = matrix(0,niter,2)
for (iter in 1:niter){
  rho.mh = rnorm(1,rho.d,sd.rho)
  sig.mh = rnorm(1,sig.d,sd.sig)
  if ((abs(rho.mh)<1) & (sig.mh>0)){
    alpha  = min(1,post(rho.mh,sig.mh)/post(rho.d,sig.d))
    if (runif(1)<alpha){
      rho.d = rho.mh
      sig.d = sig.mh
    }
  }
  draws[iter,] = c(rho.d,sig.d)
}
draws1 = draws[seq(M0+1,niter,by=L),]

par(mfrow=c(2,3),mar = c(4, 4, 1, 1))
ts.plot(draws1[,1],xlab="iterations",ylab="",main=expression(rho))
acf(draws1[,1],main="")
hist(draws1[,1],prob=TRUE,main="",xlab="")
ts.plot(draws1[,2],xlab="iterations",ylab="",main=expression(sigma))
acf(draws1[,1],main="")
hist(draws1[,2],prob=TRUE,main="",xlab="")

2.2 One-step ahead forecast

ys = matrix(0,M,n)
ys[,1] = rnorm(M,0,draws[,2])
sds = draws[,2]*sqrt((1-draws[,1]^2))
for (t in 2:n)
  ys[,t] = rnorm(M,draws[,1]*y[t-1],sds)

qys = t(apply(ys,2,quantile,c(0.05,0.5,0.95)))


par(mfrow=c(1,1))
plot(y,pch=16,ylab="Temperature anomaly (C)",ylim=range(qys[2:n,]),xlab="Time")
for (t in 2:n)
  segments(t,qys[t,1],t,qys[t,3],col=2)
points(qys[,2],pch=16,col=2,cex=0.75)

3 MODEL 2: LOCAL LEVEL MODEL

Here we model data\(=\{y_1,\ldots,y_n\}\) as a Gaussian dynamic linear model: \[\begin{eqnarray*} y_t &=& x_t + \nu_t \qquad\qquad \nu_t \sim N(0,V)\\ x_t &=& x_{t-1} + \omega_t \qquad\qquad \omega_t \sim N(0,W), \end{eqnarray*}\] where \(x_0 \sim N(m_0,C_0)\).

The prior for \((V,W)\) is \[ (V,W) \sim IG(a_v,b_v)IG(c_w,d_w), \] for \((a_v,b_v,c_w,d_w)=(2.011.01,2.01,1.01)\).

3.1 MCMC algorithm via FFBS

mcmc.joint = function(y,a1,R1,av,bv,aw,bw,M0,M,V,W){
  n      = length(y)
  V.draw = V
  W.draw = W
  niter  = M0+M
  draws  = matrix(0,niter,n+2)
  for (iter in 1:niter){
    beta.draw = ffbs(y,rep(1,n),0.0,V.draw,1.0,0.0,W.draw,a1,R1,nd=1)
    par2v     = bv+sum((y-beta.draw)^2)/2
    par2w     = bw+sum((beta.draw[2:n]-beta.draw[1:(n-1)])^2)/2
    V.draw    = 1/rgamma(1,av+n/2,par2v)
    W.draw    = 1/rgamma(1,aw+(n-1)/2,par2w)
    draws[iter,] = c(beta.draw,V.draw,W.draw)
  }
  return(draws[(M0+1):niter,])
}

ffbs = function(y,Ft,alpha,V,G,gamma,W,a1,R1,nd=1){
  n = length(y)
  if (length(Ft)==1){Ft=rep(Ft,n)}
  if (length(alpha)==1){alpha=rep(alpha,n)}
  if (length(V)==1){V=rep(V,n)}
  a = rep(0,n)
  R = rep(0,n)
  m = rep(0,n)
  C = rep(0,n)
  B = rep(0,n-1)
  H = rep(0,n-1)
  # time t=1
  a[1] = a1
  R[1] = R1
  f    = alpha[1]+Ft[1]*a[1]
  Q    = R[1]*Ft[1]**2+V[1]
  A    = R[1]*Ft[1]/Q
  m[1] = a[1]+A*(y[1]-f)
  C[1] = R[1]-Q*A**2
  # forward filtering
  for (t in 2:n){
    a[t] = gamma + G*m[t-1]
    R[t] = C[t-1]*G**2 + W
    f    = alpha[t]+Ft[t]*a[t]
    Q    = R[t]*Ft[t]**2+V[t]
    A    = R[t]*Ft[t]/Q
    m[t] = a[t]+A*(y[t]-f)
    C[t] = R[t]-Q*A**2
    B[t-1] = C[t-1]*G/R[t]
    H[t-1] = sqrt(C[t-1]-R[t]*B[t-1]**2)
  }
  # backward sampling
  theta = matrix(0,nd,n)
  theta[,n] = rnorm(nd,m[n],sqrt(C[n]))
  for (t in (n-1):1)
    theta[,t] = rnorm(nd,m[t]+B[t]*(theta[,t+1]-a[t+1]),H[t])
  if (nd==1){
    theta[1,]
  }
  else{
    theta
  }
}

# Prior for beta(0) ~ N(m0,C0)
# Prior for beta(1) ~ N(a1,R1), a1=m0,R1=C0+W
a1     = 0
R1     = 10
av     = 2.01
bv     = 1.01
aw     = 2.01
bw     = 1.01    
m0     = 0
C0     = 10

# MCMC initual values
V0     = 1
W0     = 0.5
M0     = 10000
M      = 10000

set.seed(65432)
draws = mcmc.joint(y,a1,R1,av,bv,aw,bw,M0,M,V0,W0)
qbeta = t(apply(draws[,1:n],2,quantile,c(0.05,0.5,0.95)))
draws2 = draws

3.2 MCMC-based posterior inference

par(mfrow=c(1,2))
hist(draws[,n+1],main="V",prob=TRUE,xlab="")
hist(draws[,n+2],main="W",prob=TRUE,xlab="")

par(mfrow=c(1,1))
ts.plot(qbeta,col=c(3,2,3),lwd=2)
points(y,pch=16,cex=0.5)

3.3 Comparing AR(1) and LLM

n1 = 100
par(mfrow=c(1,1))
ts.plot(qbeta[2:n1,],lwd=c(1,2,1))
lines(qys[2:n1,1],col=2,lwd=1)
lines(qys[2:n1,2],col=2,lwd=2)
lines(qys[2:n1,3],col=2,lwd=1)

4 MODEL 3: AR(1) PLUS NOISE MODEL

4.1 MCMC algorithm via FFSB

# Set up of prior hyperparameters
m0        = 0
C0        = 100
nu0       = 10
sig20     = 1
eta0      = 10
tau20     = 1
nu0sig20  = nu0*sig20
eta0tau20 = eta0*tau20
nu0n      = (nu0+n)/2
eta0n     = (eta0+n)/2

# initial values
x     = y
x0    = 0
alpha = 0
beta  = 1
# MCMC set up
thin   = 1
burnin = 1000
M      = 1000
niter  = burnin+thin*M
draws  = matrix(0,niter,n+5)
mts    = rep(0,n)
Cts    = rep(0,n)
for (iter in 1:niter){
  # Learning sig2
  sig2 = 1/rgamma(1,nu0n,(nu0sig20+sum((y-x)^2))/2)
  # Learning tau2
  xx = c(x0,x[1:(n-1)])
  tau2 = 1/rgamma(1,eta0n,(eta0tau20+sum((x-alpha-beta*xx)^2))/2)
  # Learning (alpha,beta)
  X = cbind(1,xx)
  v = solve(t(X)%*%X)
  m = v%*%t(X)%*%x
  ab = m+t(chol(v))%*%rnorm(2,0,sqrt(tau2))
  alpha = ab[1]
  beta  = ab[2]
  # Learning x0
  var  = 1/(1/C0+beta^2/tau2)
  mean = var*(m0/C0+(x[1]-alpha)*beta/tau2)
  x0   = rnorm(1,mean,sqrt(var))
  # Learning x1,...,xn
  mt = m0
  Ct = C0
  for (t in 1:n){
    at = alpha+beta*mt
    Rt = beta^2*Ct+tau2
    Qt = Rt+sig2
    At = Rt/Qt
    mt = (1-At)*at + At*y[t]
    Ct = Rt - At^2*Qt
    mts[t] = mt
    Cts[t] = Ct
  }  
  x[n] = rnorm(1,mts[n],sqrt(Cts[n]))
  for (t in (n-1):1){
    v = 1/(beta^2/tau2+1/Cts[t])
    m = v*(beta*(x[t+1]-alpha)/tau2+mts[t]/Cts[t])
    x[t] = rnorm(1,m,sqrt(v))
  }
  # Storing the draws
  draws[iter,] = c(x,sig2,tau2,alpha,beta,x0)
}
ind = seq(burnin+1,niter,by=thin)
draws3 = draws[ind,]

4.2 MCMC-based posterior inference

par(mfrow=c(4,3),mar = c(2, 2, 1, 1))
ts.plot(draws[,n+1],xlab="iterations",ylab="",main=expression(sigma^2))
acf(draws[,n+1],main="")
hist(draws[,n+1],prob=TRUE,main="",xlab="")
ts.plot(draws[,n+2],xlab="iterations",ylab="",main=expression(tau^2))
acf(draws[,n+2],main="")
hist(draws[,n+2],prob=TRUE,main="",xlab="")
ts.plot(draws[,n+3],xlab="iterations",ylab="",main=expression(alpha))
acf(draws[,n+3],main="")
hist(draws[,n+3],prob=TRUE,main="",xlab="")
ts.plot(draws[,n+4],xlab="iterations",ylab="",main=expression(beta))
acf(draws[,n+4],main="")
hist(draws[,n+4],prob=TRUE,main="",xlab="")

acf1 = acf(draws[,1],plot=FALSE)
par(mfrow=c(1,1))
plot(acf1$acf,type="l",xlab="Lag",ylab="ACF",col=grey(0.75),ylim=c(-0.2,1))
for (i in 2:n){
  acf1 = acf(draws[,i],plot=FALSE)
  lines(acf1$acf,col=grey(0.75))
}
title("ACF for each x(t)")
abline(h=0,lty=3)

4.3 Comparison

qx = t(apply(draws[,1:n],2,quantile,c(0.025,0.5,0.975)))
ts.plot(qx,xlab="Time",col=2,lty=c(2,1,2))
points(y)
lines(qbeta[,1],col=4,lty=2)
lines(qbeta[,2],col=4)
lines(qbeta[,3],col=4,lty=2)

par(mfrow=c(2,4))
hist(draws1[,1],main=expression(rho),prob=TRUE,xlab="")
hist(draws1[,2],main=expression(sigma),prob=TRUE,xlab="")
hist(draws2[,n+1],main="V",prob=TRUE,xlab="")
hist(draws2[,n+2],main="W",prob=TRUE,xlab="")
hist(draws3[,n+3],main=expression(alpha),prob=TRUE,xlab="")
hist(draws3[,n+4],main=expression(beta),prob=TRUE,xlab="")
hist(draws3[,n+1],main="V",prob=TRUE,xlab="")
hist(draws3[,n+2],main="W",prob=TRUE,xlab="")