AR(1) with Student’s \(t\) errors

Hamiltonian MC via STAN

1 Simulating some AR(1) data

We start by simulating a time-series of size \(n\) from an AR(1) process \[ y_t|y_{t-1},\mu,\phi,\sigma,\nu \sim t_\nu(\mu+\phi y_{t-1},\sigma^2), \] for \(t=1,\ldots,n\) and \(y_0=\mu/(1-\phi)\), assuming \(\phi \in (0,1)\). Below we use the values \((n,\mu,\phi,\sigma,\nu)=(200,0.1,0.9,0.1,1)\), so the errors are Cauchy.

#install.packages("rstan")
library("rstan")
options(mc.cores=4)

# Data
set.seed(3141593)
n      = 200
mu0    = 0.1
phi0   = 0.9
sigma0 = 0.1
nu0    = 1
y      = rep(mu0/(1-phi0),n)
for (t in 2:n)
  y[t] = mu0 + phi0*y[t-1] + sigma0*rt(1,df=nu0)

par(mfrow=c(1,1))
ts.plot(y)

2 Model in STAN code: AR1.stan

data {
  int<lower=0> n;   
  vector[n] y;      
}

parameters {
  real mu;                     
  real phi;  
  real<lower=0> sigma;         
  real<lower=0> nu;            
}

model {
  mu    ~ normal(0,1);
  phi   ~ normal(0.8,0.3);
  sigma ~ gamma(0.5,0.5);
  nu    ~ gamma(0.895,0.045);
  y[1]  ~ student_t(nu,0,10);
  for (t in 2:n)
    y[t] ~ student_t(nu,mu+phi*y[t-1],sigma);
}

3 Posterior inference via STAN

model = stan_model("AR1.stan")
fit   = sampling(model,list(n=n,y=y),iter=35000,warmup=10000,thin=10,chains=4)

print(fit)

params = extract(fit)
mu     = params$mu
phi    = params$phi
sigma  = params$sigma
nu     = params$nu

3.1 Checking the Markov chains

par(mfrow=c(3,4))
ts.plot(mu,xlab="Iterations",ylab="",main=expression(mu))
ts.plot(phi,xlab="Iterations",ylab="",main=expression(phi))
ts.plot(sigma,xlab="Iterations",ylab="",main=expression(sigma))
ts.plot(nu,xlab="Iterations",ylab="",main=expression(nu))
acf(mu,ylab="",main=expression(mu))
acf(phi,ylab="",main=expression(phi))
acf(sigma,ylab="",main=expression(sigma))
acf(nu,ylab="",main=expression(nu))
hist(mu,prob=TRUE,xlab="",ylab="",main=expression(mu))
abline(v=mu0,col=2,lwd=2)
hist(phi,prob=TRUE,xlab="",ylab="",main=expression(phi))
abline(v=phi0,col=2,lwd=2)
hist(sigma,prob=TRUE,xlab="",ylab="",main=expression(sigma))
abline(v=sigma0,col=2,lwd=2)
hist(nu,prob=TRUE,xlab="",ylab="",main=expression(nu))
abline(v=nu0,col=2,lwd=2)

3.2 Posterior quantiles

t(apply(cbind(mu,phi,sigma,nu),2,quantile,c(0.05,0.5,0.95)))

##               5%       50%       95%
## mu    0.07877791 0.1025356 0.1258704
## phi   0.88562498 0.9035139 0.9212185
## sigma 0.09256661 0.1139087 0.1385227
## nu    0.94907611 1.2021300 1.5250803

3.3 Bivariate marginal posteriors

par(mfrow=c(2,3))
plot(mu,phi,xlab=expression(mu),ylab=expression(phi))
plot(mu,sigma,xlab=expression(mu),ylab=expression(sigma))
plot(mu,nu,xlab=expression(mu),ylab=expression(nu))
plot(phi,sigma,xlab=expression(phi),ylab=expression(sigma))
plot(phi,nu,xlab=expression(phi),ylab=expression(nu))
plot(sigma,nu,xlab=expression(sigma),ylab=expression(nu))