# Bayesian estimation of Dynamic Regression with AR(1) parameters

I would like to draw (Bayesian) inference in a dynamic linear regression with regression parameters following independent AR(1) processes $\beta_{t,i} = \mu_i+\beta_{t-1,i}+w_{t,i}$. However, I encounter problems with my Gibbs sampler and I do not find the mistake in my approach - every comment is highly appreciated:

1. All the computations are performed using a simulated data-set. I would like to back out the AR(1) parameters of the underlying betas.

set.seed(3010)
#   Simulate autoregressive regression parameter
simAR <- function(n,sigma,mu,rho,mu0){
ts <- rep(mu0,n+1)
for (i in 1:n) ts[i+1] <- mu + rho*ts[i] + rnorm(1,sd=sqrt(sigma))
return(ts)
}
T <- 500
rhoB    <- c(0.9,0.7)
sigmaB <- rep(0.1,2)
muB <- c(-0.3,0.3)
mu0B <- rep(0,2)
beta1 <- simAR(T-1,sigmaB[1],muB[1],rhoB[1],mu0B[1])
beta2 <- simAR(T-1,sigmaB[2],muB[2],rhoB[2],mu0B[2])
beta3 <- rep(4,T)
betatrue <- cbind(beta1,beta2,beta3)
matplot(betatrue,type="l")


2. Compute the time-series $y$ with additional noise term

sigmay <- 0.3
a   <- matrix(rnorm(ncol(betatrue)*T),nc=ncol(betatrue))
y   <- matrix(apply(a*betatrue,1,sum)+rnorm(T,sd=sqrt(sigmay)))
round(rbind(OLS=coef(lm(y~a-1)),Simulated = colMeans(betatrue),True=c(muB/(1-rhoB),4)),2)
a1   a2   a3
OLS       -3.02 1.01 3.95
Simulated -2.93 1.01 4.00
True      -3.00 1.00 4.00


3. Time-varying regression with AR(1) dynamics

Investigating this model within a Bayesian regression framework is possible by setting

\begin{align*}y_t &= F_t \theta_t + v_t,& v_t \sim N(0, \sigma^2 _y) \\ \theta_t &= G \theta_{t-1} + w_t, &w_t \sim N_p (0,W) \\ \theta_0 &\sim N_p (m_0, C_0) \end{align*} where \begin{align*} F_t &=\left(X_t, 0\right) \\ \theta_t &= \left(\beta_{1,t},\beta_{2,t},\beta_{3,t},1\right) \\ G &= \left(\begin{matrix} \rho_0 & 0 & 0 & \mu_1\\ 0 & \rho_1 & 0&\mu_2\\ 0 & 0 & \rho_3 &\mu_3\\ 0&0&0&1 \end{matrix}\right)\\ W &= \left(\begin{matrix} \sigma^{2}_{1} & 0 & 0 & 0 \\ 0 & \sigma^{2} _{2} & 0 & 0\\ 0 & 0 &\sigma^{2} _{3} & 0\\ 0 & 0 & 0 & \varepsilon \end{matrix}\right)\\ \end{align*}

The Dynamic Linear Model contains the following unknown parameters $\Psi = \left(\sigma^{2} _{y} \sigma^{2} _{1}, \sigma^{2} _{2}, \sigma^{2} _{3}, \mu_1 , \mu_2 ,\mu_3 ,\rho_1 , \rho_2, \rho_3\right)$ I chose the following prior dependence structure: $\pi(\Psi) = \pi(\sigma^2 _y)\pi(W)\pi(\mu_1, \mu_2 ,\mu_3)\pi(\rho_1, \rho_2,\rho_3).$ It is convenient to define inverse-Gamma priors for the volatility terms: \begin{align*} \sigma^2_{y} & \sim \mathcal{G}^{-1}\left(c_y,C_y\right) \\ \sigma^2_{i} & \sim \mathcal{G}^{-1}\left(c_i,C_i\right) \end{align*} The conditional priors for the drift coefficients are chosen to be normal $\mu_i |\sigma^2 _i \sim N\left(b_0,\frac{\sigma^2_{i}}{N_0}\right).$ The AR(1) parameters ($\rho_i$) are not bounded to the stationary region by theprior choice as I set $\pi(\rho_i) \sim N(\psi_0,\tau_0) = N(0, 1)$.

4. Sampling procedure

The Gibbs sampling consists of the following steps:

• Start with initial guess $\Psi^0$
• Sample the states $\Theta^i$ from $\pi\left(\Theta|\Psi^i,Y\right)$
• Generate draws $\Psi^{i+1}$ from $\pi\left(\Psi|\Theta^i ,Y\right)$

Sampling the states is done via Forward Filtering Backward Sampling. Conditional on having sampled $\Theta^i$, the unknown variances $\sigma_{y}^{2}$ and $\sigma_{i}^{2}$ are computed using $\sigma^2_{y}| \dots \sim \mathcal{G}^{-1}\left(c_y + \frac{T}{2}, C_y + \frac{1}{2}\sum\limits_{t=1}^T (y_t-(F_t \theta_t))^2\right)$ and $\sigma^2_i| \dots \sim \mathcal{G}^{-1}\left(c_i + \frac{T+1}{2},\frac{1}{2}\sum\limits_{t=1}^T (\theta_{i,t}-(G\theta_{t-1})_i+ \frac{N_0}{2}\left(b_0 - \mu_i\right)^2+ C_i\right).$

For the unknown drifts $\mu_i$ I obtain $p(\mu_i|\dots) \sim N\left(\frac{1}{T+N_0}\left(\sum\limits_{t=1}^T \left(\theta_{t,i}-\rho_i\theta_{t-1,i}\right)+N_0 b_0\right), \sigma^2 _i /(N_0 + T)\right)$.

The AR parameters $\rho_i$ are computed using

$\rho_i | \dots\sim N\left(\psi_j^T,\tau_j ^T\right)$ with $\tau_j ^T := \left(\frac{1}{\sigma_i ^2}\sum\limits_{t=1}^T {\theta_{t-1} ^i} ^2 +\frac{1}{\tau^2}\right)^{-1}$ and $\psi_j^T := \left(\tau_j ^T\right)^{-1} \left[\frac{1}{\sigma_i ^2 }\sum\limits_{t=1}^T \left(\theta_t ^i - \mu_i\right)\theta_{t-1} ^i + \frac{1}{\tau ^2} \psi\right]$.

5. Simulation The R Code (using the DLM package) is constructed as follows (and in large parts borrowed from this excellent Book (code provided on home page).

library(dlm)

p       <- 4            # State dimension
m       <- ncol(y)
T       <- nrow(y)
B       <- function(phi.init,Mu.init,tot=p) diag(c(phi.init,rep(1,p-length(phi.init))))+ cbind(matrix(0,nc=p-1,nr=p),c(Mu.init,0))

a2       <- cbind(a,0)
psi0       <- 0
tau0       <- 1
shapeY     <- 2
rateY      <- 0.1
shapeTheta <- 2
rateTheta  <- 0.1
Mu0          <- 0
N0       <- 10

######MCMC - Initialization############################

burn        <- 10000
MC      <- 10000 + burn
gibbsTheta  <- array(NA, dim=c(MC,T+1,p))
gibbsPsi    <- matrix(NA, nrow=MC, ncol=p-1)
gibbsV  <- matrix(NA, nrow=MC, ncol=m)
gibbsW  <- matrix(NA, nrow=MC, ncol=p-1)
gibbsMu <- matrix(NA, nrow=MC, ncol=p-1)

tiny        <- 1e-12

phi.init    <- rnorm(p-1,mean=psi0,sd=tau0)
V.init  <- 1/rgamma(1,shapeY,rateY)
W.init  <- 1/rgamma(1,shapeTheta,rateTheta)
Mu.init     <- rnorm(p-1,mean=0,sd=W.init/N0)
mod         <- dlm(m0=c(rep(0,p-1),rep(1,1)),
C0=diag(c(rep(1e01,p-1),tiny)),
FF=matrix(1,nrow=1,ncol=p),
JFF=matrix(1,nc=p,nr=1),
V= var(residuals(lm(y~a2-1))),
GG=B(phi.init,Mu.init),
W=diag(c(rep(W.init,p-1),tiny)),X=a2)

######Gibbs Sampling###################################

for (i in 1:MC){
if(i%%100==0){ cat('# Iteration', i,'- V:     ',round(median(gibbsV[1:i],na.rm=TRUE),3),'\n'); flush.console()}
# generate the states by FFBS
modFilt         <- dlmFilter(y, mod, simplify=TRUE)
theta       <- dlmBSample(modFilt)
gibbsTheta[i,,]     <- theta
# generate the W_j
theta.center    <- theta[-1,] - t(mod$GG %*% t(theta[-(T+1),])) SStheta <- apply((theta.center)^2, 2,sum) phiTheta <- rep(NA,p-1) for(j in 1:(p-1)){phiTheta[j]<-rgamma(1, shape=shapeTheta+(T+1)/2, rate=rateTheta+SStheta[j]/2+N0/2*(Mu.init[j]-Mu0)^2)} gibbsW[i,] <- 1 / phiTheta mod$W        <- diag(c(gibbsW[i,],tiny))
# generate the Mu_j
MuMean <- rep(NA,p-1)
for(j in 1:(p-1)) MuMean[j] <- 1/(T+N0)*(N0*Mu0+sum(theta[-1,j]-mod$GG[j,j]*theta[-T,j])) MuVar <- gibbsW[i,]/(T+N0) for(j in 1:(p-1)) {Mu.init[j] <- rnorm(1,mean=MuMean[j],sd=sqrt(MuVar[j]))} gibbsMu[i,] <- Mu.init mod$GG     <- B(phi.init,gibbsMu[i,])
# generate the V_i
y.center    <- rep(NA,T)
for(t in 1:T) y.center[t] <- y[t] - mod$X[t,]%*%theta[t+1,] SSy <- sum(y.center^2) gibbsV[i,] <- 1/rgamma(m, shape=shapeY+T/2, rate=rateY+SSy/2) mod$V  <- gibbsV[i,]
# generate the AR parameters psi_1, psi_2
phi.init    <- rep(NA,p-1)
for (j in 1:(p-1)){
tau <- 1/(1/tau0 + 1/gibbsW[i,j] * crossprod(theta[-1,j]))
psi <- tau * (1/gibbsW[i,j] * sum((theta[-1,j]-gibbsMu[i,j])*theta[-T,j]) + psi0/tau0)
phi.init[j] <- rnorm(1, psi, sd=sqrt(tau))
}
gibbsPsi[i,] <- phi.init
mod\$GG       <- B(gibbsPsi[i,],gibbsMu[i,])
}


6. Cleaning, Thining and evaluation After the simulation steps I obtain the following output for the states (I obtain every 5th observation to reduce autocorrelation in the sample).

 # Remove burn-in phase + Thin-out
thinfactor <- 5
thin <- function(dat,fact=thinfactor) as.matrix(dat)[seq(from=1,to=nrow(as.matrix(dat)),by=fact),]

gibbsW  <- thin(gibbsW[-(1:burn),])
gibbsMu     <- thin(gibbsMu[-(1:burn),])
gibbsPsi    <- thin(gibbsPsi[-(1:burn),])
gibbsTheta  <- gibbsTheta[seq(from=burn+1,to=MC-burn,by=thinfactor),,]
gibbsV  <- thin(gibbsV[-(1:burn),])

# Theta (state)
postmeanTheta   <- apply(gibbsTheta,c(2,3),function(x) mean(x,na.rm=TRUE))
postquantileTheta <- apply(gibbsTheta,c(2,3),function(x)quantile(x,c(0.05,0.95),na.rm=TRUE))
par(mar = c(2, 4, 1, 1) + 0.1, cex = 0.6, mfrow=c(1,3))

for(j in 1:3){
ytmp   <- postmeanTheta[-(1:5),j]
y01 <- postquantileTheta[1,-(1:5),j]
y09 <- postquantileTheta[2,-(1:5),j]
plot(ytmp,type="l",ylim=range(ytmp,y01,y09),main=paste('Posterior distribution of Theta',j))
lines(y01,col="grey")
lines(y09,col="grey")
lines(betatrue[,j],col='red')
}


Apparently, the posterior distribution is not really capturing what I aimed at. The plot shows the posterior draws of the states and should (at least this is what I expect) be somewhat comparable to the 'true' states which are plotted in red. Why isn't this the case?

• Are there too much parameters? (I doubt so)
• Is there an error in my Gibbs sampler derivations (I can provide the full calculations if needed, but the computations do not seem odd to me)
• Is there an implementation issue?
• Is this whole simulation study not suited to test whether my code works?