3
$\begingroup$

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")

Time-series of 'true' betas

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')
    }

enter image description here

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?
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.