Gibbs Sampler: Dynamic CAR Model

Author

Affiliation

Elliot Shannon

Michigan State University

Setting

We consider data observed over a spatial domain which is partitioned into a fixed number of areal units. The data are observed for each areal unit across discrete time steps. We aim to model these data using spatially and temporally explicit covariates, as well as a spatially-varying intercept term which is allowed to evolve dynamically over time. We will employ conditional autoregressive (CAR) spatial structures for these random effects. Examples of data such as described here include annual disease prevalence for US counties or quarterly average home values for census tracts in Michigan.

Example of a CAR spatial effect evolving over time (from left to right).

Model

For areal units \(j = 1, \ldots, J\) and time points \(t = 1, \ldots, T\), let

\(y_{j,t} \in \mathbb{R}\) be the response value.
\(\textbf{x}_{j, t} = ( 1, x_{1,j,t} \cdots, x_{P - 1,j,t}) \in \mathbb{R}^P\) be the vector of covariates (predictor variables).
\(\boldsymbol{\beta}= \begin{pmatrix} \beta_{1} \\ \vdots \\ \beta_{P} \end{pmatrix} \in \mathbb{R}^P\) be the vector of regression coefficients.
\(u_{j, t} \in \mathbb{R}\) be the spatio-temporally varying random intercept.
\(\varepsilon_{j, t} \in \mathbb{R}\) be the residual error term.

For areal unit \(j\) at time \(t\), the proposed model/priors is then

\[\begin{equation}\label{eq:mod} \begin{aligned} y_{j,t} &= u_{j,t} + \textbf{x}_{j,t} \boldsymbol{\beta}+ \varepsilon_{j,t}, \quad \varepsilon_{j,t} \sim N(0, \sigma^2), \\ u_{j, t} &= u_{j, t-1} + w_{j, t}, \quad u_{\cdot, 0} \equiv 0 \\ \sigma^2 &\sim IG(a_\sigma, b_\sigma), \quad \boldsymbol{\beta}\sim MVN( \boldsymbol{\mu} _\beta, \boldsymbol{\Sigma} _\beta) \end{aligned} \end{equation}\]

where \(\textbf{w}_{t} = (w_{1, t}, \ldots, w_{J, t})^\top\) follows from a conditional autoregressive (CAR) spatial structure as

\[\begin{equation} \begin{aligned}\label{eq:car} \textbf{w}_t &\sim MVN(0, \tau^2_t \textbf{Q}(\rho_t)), \quad \textbf{Q}(\rho_t) = (\textbf{D}- \rho_t \textbf{W})^{-1} \\ \tau^2_t &\sim IG(a_t, b_t), \quad \rho_t \sim U(c_t, d_t). \end{aligned} \end{equation}\]

Here, \(\tau^2_t\) is the spatial variance parameter, \(\textbf{Q}(\rho_t)\) is the CAR correlation matrix depending on spatial dependence parameter \(\rho_t\), \(\textbf{W}\) is a \(J \times J\) binary spatial adjacency matrix, and \(\textbf{D}\) is a \(J \times J\) diagonal matrix whose elements correspond to the row sums of \(\textbf{W}\). We assume \(\textbf{W}\) and \(\textbf{D}\) to be fixed over time, but allow for unique variance and spatial dependence parameters for each time step.

Stacking

Stacking these terms for each time \(t = 1, \ldots, T\), we can also write them in the form

\(\textbf{y}_{t} = \begin{pmatrix} y_{1, t} \\ \vdots \\ y_{J,t} \end{pmatrix} \in \mathbb{R}^J\)
\(\textbf{X}_{t} = \begin{pmatrix} \textbf{x}_{1,t} \\ \vdots \\ \textbf{x}_{J,t} \end{pmatrix} = \begin{pmatrix} 1 & x_{1,1,t} & \cdots & x_{P - 1,1,t} \\ & & \vdots \\ 1 & x_{1,J,t} & \cdots & x_{P - 1,J,t} \end{pmatrix} \in \mathbb{R}^{J \times P}\)
\(\textbf{u}_{t} = \begin{pmatrix} u_{1, t} \\ \vdots \\ u_{J,t} \end{pmatrix} \in \mathbb{R}^J\)
\( \boldsymbol{\varepsilon} _{t} = \begin{pmatrix} \varepsilon_{1, t} \\ \vdots \\ \varepsilon_{J,t} \end{pmatrix} \in \mathbb{R}^J\)

and the model for time \(t\) is then written as

\[\begin{equation}\label{eq:mod_stacked} \begin{aligned} \textbf{y}_{t} &= \textbf{u}_{t} + \textbf{X}_{t} \boldsymbol{\beta}+ \boldsymbol{\varepsilon} _{t}, \quad \boldsymbol{\varepsilon} _{t} \sim MVN(\textbf{0}, \sigma^2 \textbf{I}_J), \\ \textbf{u}_{t} &= \textbf{u}_{t-1} + \textbf{w}_{t}, \quad \textbf{u}_{0} \equiv 0 \\ \sigma^2 &\sim IG(a_\sigma, b_\sigma), \quad \boldsymbol{\beta}\sim MVN( \boldsymbol{\mu} _\beta, \boldsymbol{\Sigma} _\beta) \\ \textbf{w}_t &\sim MVN(0, \tau^2_t \textbf{Q}(\rho_t)), \quad \textbf{Q}(\rho_t) = (\textbf{D}- \rho_t \textbf{W})^{-1} \\ \tau^2_t &\sim IG(a_t, b_t), \quad \rho_t \sim U(c_t, d_t). \end{aligned} \end{equation}\]

where \(\textbf{I}_J\) is a \(J \times J\) identity matrix with diagonal elements equal to \(1\), and \(\textbf{0}\) is a vector of \(0\)’s.

Here, \(N(\mu, \sigma^2)\) is the normal distribution with mean \(\mu \in \mathbb{R}\), variance \(\sigma^2 > 0\), and pdf

\[\begin{equation}\label{eq:normal} p(x \mid \mu, \sigma^2) = (2 \pi \sigma^2)^{-1/2} \exp\left(\frac{-1}{2\sigma^2} (x - \mu)^2 \right) \end{equation}\]

\(MVN( \boldsymbol{\mu} , \boldsymbol{\Sigma} )\) is the multivariate normal distribution with mean vector \( \boldsymbol{\mu} \in \mathbb{R}^N\), covariance matrix \( \boldsymbol{\Sigma} \in \mathbb{R}^{N \times N}\), and pdf

\[\begin{equation}\label{eq:mvn} p(\textbf{x}\mid \boldsymbol{\mu} , \boldsymbol{\Sigma} ) = (2 \pi)^{-N/2} | \boldsymbol{\Sigma} |^{-1/2} \exp\left(- \frac{1}{2} (\textbf{x}- \boldsymbol{\mu} )^\top \boldsymbol{\Sigma} ^{-1} (\textbf{x}- \boldsymbol{\mu} )\right), \end{equation}\]

\(IG(a, b)\) is the Inverse Gamma distribution with shape parameter \(a > 0\), scale parameter \(b > 0\), and pdf

\[\begin{equation}\label{eq:ig} p(\sigma^2 \mid a, b) = \frac{b^a}{\Gamma (a)} (\sigma^2)^{-a - 1} \exp\left( \frac{-b}{\sigma^2} \right), \end{equation}\]

and \(U(c, d)\) is the continuous Uniform distribution with support \([c, d]\) and pdf

\[\begin{equation}\label{eq:unif} p(\rho \mid c, d) = \begin{cases} \frac{1}{d - c} & \text{if } c \leq \rho \leq d \\ 0 & \text{otherwise}. \end{cases} \end{equation}\]

Given this completed model specification, we would like to conduct inference on the model parameters given data \(\textbf{y}_t\) and \(\textbf{X}_t\), \(t = 1, \ldots, T\). These parameters include \(\boldsymbol{\beta}, \textbf{u}_1, \ldots, \textbf{u}_T, \tau^2_1, \ldots, \tau^2_T, \rho_1, \ldots, \rho_T, \sigma^2\). We assume that prior parameters are given (specified by the user).

Posterior

We know that the posterior distribution in Bayesian inference is proportional to the likelihood (data distribution) times the prior. Here, our likelihood is the distribution of the data \(\textbf{y}_t, t = 1, \ldots, T\). Since the distribution of the residuals is \(MVN\), we will have a \(MVN\) likelihood, written as

\[\begin{equation}\label{eq:likelihood} \prod_{t = 1}^T MVN(\textbf{y}_t \mid \textbf{u}_t + \textbf{X}_t \boldsymbol{\beta}, \sigma^2 \textbf{I}_J). \end{equation}\]

We will also have priors (some of which are induced) for each parameter for which we seek inference. Specifically, we have prior distributions of the form

\[\begin{equation}\label{eq:priors} \begin{aligned} &\prod_{t = 1}^T MVN(\textbf{u}_t \mid \textbf{u}_{t - 1}, \tau^2_t \textbf{Q}(\rho_t)) \times MVN(\boldsymbol{\beta}\mid \boldsymbol{\mu} _\beta, \boldsymbol{\Sigma} _\beta) \times IG(\sigma^2 \mid a, b) \\ &\quad \times \prod_{t = 1}^T IG(\tau^2_t \mid a_t, b_t) \times \prod_{t = 1}^T U(\rho_t \mid c_t, d_t) \end{aligned} \end{equation}\]

Given (\(\ref{eq:likelihood}\)) and (\(\ref{eq:priors}\)), their product will be proportional to the posterior distribution, and is given as

\[\begin{equation}\label{eq:posterior} \begin{aligned} &\prod_{t = 1}^T MVN(\textbf{y}_t \mid \textbf{u}_t + \textbf{X}_t \boldsymbol{\beta}, \sigma^2 \textbf{I}_J) \times \prod_{t = 1}^T MVN(\textbf{u}_t \mid \textbf{u}_{t - 1}, \tau^2_t \textbf{Q}(\rho_t)) \times MVN(\boldsymbol{\beta}\mid \boldsymbol{\mu} _\beta, \boldsymbol{\Sigma} _\beta) \\ &\quad \times IG(\sigma^2 \mid a, b) \times \prod_{t = 1}^T IG(\tau^2_t \mid a_t, b_t) \times \prod_{t = 1}^T U(\rho_t \mid c_t, d_t) \end{aligned} \end{equation}\]

Moving forward, we only work with this proportional posterior (\(\ref{eq:posterior}\)).

Working it out

We will now go through each of the parameters in model (\(\ref{eq:mod}\)) and derive their full conditional distributions. Moving forward, we will use different colors to help keep track of different prior and likelihood terms as we work through the derivations. This will help us see where each component of the full conditional distribution is coming from.

Update \(\boldsymbol{\beta}\)

Looking at (\(\ref{eq:posterior}\)), we see that \(\boldsymbol{\beta}\) appears in \(T + 1\) terms, namely in and in its . The product of these terms looks like

\[\begin{equation*} \textcolor{Emerald}{\prod_{t = 1}^T MVN(\textbf{y}_t \mid \textbf{u}_t + \textbf{X}_t \boldsymbol{\beta}, \sigma^2 \textbf{I}_J)} \times \textcolor{OrangeRed}{MVN(\boldsymbol{\beta}\mid \boldsymbol{\mu} _\beta, \boldsymbol{\Sigma} _\beta)}. \end{equation*}\]

This is the product of multivariate normal distributions, so we know through conjugacy that the full conditional distribution will also be multivariate normal. Therefore, we simply need to identify the mean and covariance matrix of this full conditional distribution, which can be obtained by .

Incompleting the square

See section 4.2 in James Clark’s book ‘Models for Ecological Data’ for more details.

Following (\(\ref{eq:mvn}\)), we know that the kernel of the \(MVN\) distribution is of the form

\[\begin{equation*} p(\textbf{x}\mid \boldsymbol{\mu} , \boldsymbol{\Sigma} ) \propto \exp\left(- \frac{1}{2} (\textbf{x}- \boldsymbol{\mu} )^\top \boldsymbol{\Sigma} ^{-1} (\textbf{x}- \boldsymbol{\mu} )\right) \end{equation*}\]

where \(\propto\) mean “proportional to”. This structure gives us some insight into the placement of the parameters in the pdf. To see this further, we can focus on the terms inside the exponential, ignoring the \(-\frac{1}{2}\) coefficient. Expanding this term gives us:

\[\begin{align}\label{eq:expand_mvn} (\textbf{x}- \boldsymbol{\mu} )^\top \boldsymbol{\Sigma} ^{-1} (\textbf{x}- \boldsymbol{\mu} ) &= \textbf{x}^\top \boldsymbol{\Sigma} ^{-1} \textbf{x}- \textbf{x}^\top \boldsymbol{\Sigma} ^{-1} \boldsymbol{\mu} - \boldsymbol{\mu} ^\top \boldsymbol{\Sigma} ^{-1} \textbf{x}+ \boldsymbol{\mu} ^\top \boldsymbol{\Sigma} ^{-1} \boldsymbol{\mu} \end{align}\]

Now comes the trick. Writing the mean as a function of the covariance matrix, we can define

\[\begin{align*} \boldsymbol{\mu} &= \textbf{V}\textbf{v}\\ \boldsymbol{\Sigma} &= \textbf{V} \end{align*}\]

This allows us to rewrite (\(\ref{eq:expand_mvn}\)) as

\[\begin{align*} \textbf{x}^\top \textbf{V}^{-1} \textbf{x}&- \textbf{x}^\top \textbf{V}^{-1} \textbf{V}\textbf{v}- \textbf{v}^\top \textbf{V}^\top \textbf{V}^{-1} \textbf{x}+ \textbf{v}^\top \textbf{V}^\top \textbf{V}^{-1} \textbf{V}\textbf{v}\\ &= \textbf{x}^\top \textbf{V}^{-1} \textbf{x}- \textbf{x}^\top \textbf{v}- \textbf{v}^\top \textbf{x}+ \textbf{v}^\top \textbf{V}\textbf{v}\\ &= \textbf{x}^\top \textbf{V}^{-1} \textbf{x}- \textbf{x}^\top \textbf{v}\cdots \end{align*}\]

This is super helpful, because when we are mashing together normal priors and likelihoods, we can combine like terms and look for \(\textbf{V}^{-1}\) and \(\textbf{v}\) to get the mean and variance of the full conditional \(MVN\). They will always be in the same place! \(\textbf{V}^{-1}\) will be sandwiched between the variable of interest (i.e. the \(\textbf{x}^\top \textbf{V}^{-1} \textbf{x}\) term), and \(\textbf{v}\) will be attached to the negative transpose of the variable of interest (i.e. \(- \textbf{x}^\top \textbf{v}\)) (in this example, \(\textbf{x}\) is the variable of interest). Once we identify \(\textbf{V}^{-1}\) and \(\textbf{v}\), we can sample from the full conditional distribution as \(MVN(\textbf{V}\textbf{v}, \textbf{V})\).

Looking back to the full conditional for \(\boldsymbol{\beta}\), we will find \(\textbf{V}\) and \(\textbf{v}\) in the following terms of the product

\[\begin{equation*} \boldsymbol{\beta}^\top \textbf{V}^{-1} \boldsymbol{\beta}- \boldsymbol{\beta}^\top \textbf{v}\cdots. \end{equation*}\]

So, all we have to do is multiply \(T + 1\) multivariate normal pdfs (see equation \(\ref{eq:mvn}\)), combine like terms in \(\boldsymbol{\beta}\), and identify \(\textbf{V}^{-1}\) and \(\textbf{v}\) from where they appear above. Since we are incompleting the square, we don’t actually need to multiply the complete pdfs, just the exponential parts where \(\textbf{V}\) and \(\textbf{v}\) will appear. Ignoring the other components of the pdfs, we have

\[\begin{align*} &\textcolor{Emerald}{\prod_{t = 1}^T MVN(\textbf{y}_t \mid \textbf{u}_t + \textbf{X}_t \boldsymbol{\beta}, \sigma^2 \textbf{I}_J)} \times \textcolor{OrangeRed}{MVN(\boldsymbol{\beta}\mid \boldsymbol{\mu} _\beta, \boldsymbol{\Sigma} _\beta)} \\ &\propto \textcolor{Emerald}{\prod_{t = 1}^T \exp \left( (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta})^\top (\sigma^2 \textbf{I}_J)^{-1} (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta})\right)} \\ &\quad \times \textcolor{OrangeRed}{\exp \left( (\boldsymbol{\beta}- \boldsymbol{\mu} _\beta)^\top \boldsymbol{\Sigma} _\beta^{-1} (\boldsymbol{\beta}- \boldsymbol{\mu} _\beta)\right)} \\ &= \exp \left( \textcolor{Emerald}{\sum_{t=1}^T \left( (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta})^\top (\sigma^2 \textbf{I}_J)^{-1} (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta}) \right)} + \textcolor{OrangeRed}{(\boldsymbol{\beta}- \boldsymbol{\mu} _\beta)^\top \boldsymbol{\Sigma} _\beta^{-1} (\boldsymbol{\beta}- \boldsymbol{\mu} _\beta)} \right) \\ &= \exp \Bigg( \textcolor{Emerald}{\sum_{t=1}^T \left( \textbf{y}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{y}_t - \textbf{y}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{u}_t - \textbf{y}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{X}_t \boldsymbol{\beta}\right.} \big. \\ &\quad \quad \quad \quad \left. \textcolor{Emerald}{ - \textbf{u}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{y}_t + \textbf{u}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{u}_t + \textbf{u}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{X}_t \boldsymbol{\beta}} \right. \\ &\quad \quad \quad \quad \textcolor{Emerald}{\left. - \boldsymbol{\beta}^\top \textbf{X}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{y}_t + \boldsymbol{\beta}^\top \textbf{X}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{u}_t + \boldsymbol{\beta}^\top \textbf{X}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{X}_t \boldsymbol{\beta}\right)} \\ &\quad \quad \quad \quad \big. \textcolor{OrangeRed}{+ \boldsymbol{\beta}^\top \boldsymbol{\Sigma} _\beta^{-1} \boldsymbol{\beta}- \boldsymbol{\beta}^\top \boldsymbol{\Sigma} _\beta^{-1} \boldsymbol{\mu} _\beta - \boldsymbol{\mu} _\beta \boldsymbol{\Sigma} _\beta^{-1} \boldsymbol{\beta}+ \boldsymbol{\mu} _\beta^\top \boldsymbol{\Sigma} _\beta^{-1} \boldsymbol{\mu} _\beta} \Bigg) \\ &\propto \exp \left( \boldsymbol{\beta}^\top \left( \textcolor{Emerald}{\sum_{t=1}^T \textbf{X}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{X}_t} + \textcolor{OrangeRed}{ \boldsymbol{\Sigma} _\beta^{-1}} \right) \boldsymbol{\beta}- \boldsymbol{\beta}^\top \left( \textcolor{Emerald}{\sum_{t=1}^T \textbf{X}_t^\top (\sigma^2 \textbf{I}_J)^{-1} (\textbf{y}_t - \textbf{u}_t)} + \textcolor{OrangeRed}{ \boldsymbol{\Sigma} _\beta^{-1} \boldsymbol{\mu} _\beta} \right) - \cdots \right) \end{align*}\]

So we can identify

\[\begin{align*} &\textbf{V}^{-1} = \textcolor{Emerald}{\sum_{t=1}^T \textbf{X}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{X}_t} + \textcolor{OrangeRed}{ \boldsymbol{\Sigma} _\beta^{-1}} \\ &\textbf{v}= \textcolor{Emerald}{\sum_{t=1}^T \textbf{X}_t^\top (\sigma^2 \textbf{I}_J)^{-1} (\textbf{y}_t - \textbf{u}_t)} + \textcolor{OrangeRed}{ \boldsymbol{\Sigma} _\beta^{-1} \boldsymbol{\mu} _\beta} \end{align*}\]

and we can update \(\boldsymbol{\beta}\) from its full conditional as \(MVN(\textbf{V}\textbf{v}, \textbf{V})\).

Update \(\textbf{u}_t, t = 1, \ldots, T-1\)

Looking at (\(\ref{eq:posterior}\)), we see that the term \(\textbf{u}_t\) appears in \(3\) terms, namely the , the , and in the . The product of these terms is the product of multivariate normal distributions,

\[\begin{equation*} \textcolor{Emerald}{MVN(\textbf{y}_t \mid \textbf{u}_t + \textbf{X}_t \boldsymbol{\beta}, \sigma^2 \textbf{I}_J)} \times \textcolor{RoyalBlue}{MVN(\textbf{u}_t \mid \textbf{u}_{t-1}, \tau^2_t \textbf{Q}(\rho_t))} \times \textcolor{OliveGreen}{MVN(\textbf{u}_{t+1} \mid \textbf{u}_{t}, \tau^2_{t+1} \textbf{Q}(\rho_{t+1}))}. \end{equation*}\]

Again, we will to identify the mean and covariance matrix of the resulting multivariate normal full conditional distribution for \(\textbf{u}_t\). Specifically, we will find \(\textbf{V}\) and \(\textbf{v}\) in the following terms of the product

\[\begin{equation*} \textbf{u}_t^\top \textbf{V}^{-1} \textbf{u}_t - \textbf{u}_t^\top \textbf{v}\cdots. \end{equation*}\]

We then have

\[\begin{align*} &\textcolor{Emerald}{MVN(\textbf{y}_t \mid \textbf{u}_t + \textbf{X}_t \boldsymbol{\beta}, \sigma^2 \textbf{I}_J)} \times \textcolor{RoyalBlue}{MVN(\textbf{u}_t \mid \textbf{u}_{t-1}, \tau^2_t \textbf{Q}(\rho_t))} \times \textcolor{OliveGreen}{MVN(\textbf{u}_{t+1} \mid \textbf{u}_{t}, \tau^2_{t+1} \textbf{Q}(\rho_{t+1}))} \\ &\propto \textcolor{Emerald}{\exp \left( (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta})^\top (\sigma^2 \textbf{I}_J)^{-1} (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta}) \right)} \times \textcolor{RoyalBlue}{\exp \left( (\textbf{u}_t - \textbf{u}_{t-1})^\top (\tau^2_t \textbf{Q}(\rho_t))^{-1} (\textbf{u}_t - \textbf{u}_{t-1}) \right)} \\ &\quad \times \textcolor{OliveGreen}{\exp \left( (\textbf{u}_{t+1} - \textbf{u}_{t})^\top (\tau^2_{t+1} \textbf{Q}(\rho_{t+1}))^{-1} (\textbf{u}_{t+1} - \textbf{u}_{t}) \right)} \\ &= \exp \left( \textcolor{Emerald}{(\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta})^\top (\sigma^2 \textbf{I}_J)^{-1} (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta})} + \textcolor{RoyalBlue}{(\textbf{u}_t - \textbf{u}_{t-1})^\top (\tau^2_t \textbf{Q}(\rho_t))^{-1} (\textbf{u}_t - \textbf{u}_{t-1})} \right. \\ &\quad \quad \quad \quad + \left. \textcolor{OliveGreen}{(\textbf{u}_{t+1} - \textbf{u}_{t})^\top (\tau^2_{t+1} \textbf{Q}(\rho_{t+1}))^{-1} (\textbf{u}_{t+1} - \textbf{u}_{t})} \right) \\ &= \exp \left( \textcolor{Emerald}{\textbf{y}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{y}_t - \textbf{y}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{u}_t - \textbf{y}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{X}_t \boldsymbol{\beta}} \right. \\ &\quad \quad \quad \quad \textcolor{Emerald}{- \textbf{u}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{y}_t + \textbf{u}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{u}_t + \textbf{u}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{X}_t \boldsymbol{\beta}} \\ &\quad \quad \quad \quad \textcolor{Emerald}{- \boldsymbol{\beta}^\top \textbf{X}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{y}_t + \boldsymbol{\beta}^\top \textbf{X}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{u}_t + \boldsymbol{\beta}^\top \textbf{X}_t^\top (\sigma^2 \textbf{I}_J)^{-1} \textbf{X}_t \boldsymbol{\beta}} \\ &\quad \quad \quad \quad \textcolor{RoyalBlue}{+ \textbf{u}_t^\top (\tau_t^2 \textbf{Q}(\rho_t))^{-1} \textbf{u}_t - \textbf{u}_t^\top (\tau_t^2 \textbf{Q}(\rho_t))^{-1} \textbf{u}_{t-1} - \textbf{u}_{t-1}^\top (\tau_t^2 \textbf{Q}(\rho_t))^{-1} \textbf{u}_t + \textbf{u}_{t-1}^\top (\tau_t^2 \textbf{Q}(\rho_t))^{-1} \textbf{u}_{t-1} } \\ &\quad \quad \quad \quad \textcolor{OliveGreen}{+ \textbf{u}_{t+1}^\top (\tau_{t+1}^2 \textbf{Q}(\rho_{t+1}))^{-1} \textbf{u}_{t+1} - \textbf{u}_{t+1}^\top (\tau_{t+1}^2 \textbf{Q}(\rho_{t+1}))^{-1} \textbf{u}_{t}} \\ &\quad \quad \quad \quad \left. \textcolor{OliveGreen}{- \textbf{u}_{t}^\top (\tau_{t+1}^2 \textbf{Q}(\rho_{t+1}))^{-1} \textbf{u}_{t+1} + \textbf{u}_{t}^\top (\tau_{t+1}^2 \textbf{Q}(\rho_{t+1}))^{-1} \textbf{u}_{t}} \right) \\ &\propto \exp \left( \textbf{u}_t^\top \left( \textcolor{Emerald}{(\sigma^2 \textbf{I}_J)^{-1}} + \textcolor{RoyalBlue}{(\tau_t^2 \textbf{Q}(\rho_t))^{-1}} + \textcolor{OliveGreen}{(\tau_{t+1}^2 \textbf{Q}(\rho_{t+1}))^{-1}} \right) \textbf{u}_t \right. \\ &\quad \quad \left. - \textbf{u}_t^\top \left( \textcolor{Emerald}{(\sigma^2 \textbf{I}_J)^{-1}(\textbf{y}_t - \textbf{X}_t \boldsymbol{\beta})} + \textcolor{RoyalBlue}{(\tau_t^2 \textbf{Q}(\rho_t))^{-1}\textbf{u}_{t-1}} + \textcolor{OliveGreen}{(\tau_{t+1}^2 \textbf{Q}(\rho_{t+1}))^{-1}\textbf{u}_{t+1}} \right) - \cdots \right) \end{align*}\]

So we can identify

\[\begin{align*} \textbf{V}^{-1} &= \textcolor{Emerald}{(\sigma^2 \textbf{I}_J)^{-1}} + \textcolor{RoyalBlue}{(\tau_t^2 \textbf{Q}(\rho_t))^{-1}} + \textcolor{OliveGreen}{(\tau_{t+1}^2 \textbf{Q}(\rho_{t+1}))^{-1}} \\ \textbf{v}&= \textcolor{Emerald}{(\sigma^2 \textbf{I}_J)^{-1}(\textbf{y}_t - \textbf{X}_t \boldsymbol{\beta})} + \textcolor{RoyalBlue}{(\tau_t^2 \textbf{Q}(\rho_t))^{-1}\textbf{u}_{t-1}} + \textcolor{OliveGreen}{(\tau_{t+1}^2 \textbf{Q}(\rho_{t+1}))^{-1}\textbf{u}_{t+1}} \end{align*}\]

and we can update \(\textbf{u}_t\) from its full conditional as \(MVN(\textbf{V}\textbf{v}, \textbf{V})\).

Update \(\textbf{u}_T\)

To derive the full conditional for \(\textbf{u}_T\), it will be the same as above, except the terms for will not be present, and the full conditional mean and covariance terms are of the form

\[\begin{align*} \textbf{V}^{-1} &= \textcolor{Emerald}{(\sigma^2 \textbf{I}_J)^{-1}} + \textcolor{RoyalBlue}{(\tau_T^2 \textbf{Q}(\rho_T))^{-1}} \\ \textbf{v}&= \textcolor{Emerald}{(\sigma^2 \textbf{I}_J)^{-1}(\textbf{y}_T - \textbf{X}_T \boldsymbol{\beta})} + \textcolor{RoyalBlue}{(\tau_T^2 \textbf{Q}(\rho_T))^{-1}\textbf{u}_{T-1}} \end{align*}\]

and we can update \(\textbf{u}_T\) from its full conditional as \(MVN(\textbf{V}\textbf{v}, \textbf{V})\).

Alternative update for \(u_{j,t}, j = 1, \ldots, J, \quad t = 1, \ldots, T\)

Previously, we have jointly updated \(\textbf{u}_t\) as a multivariate normal distribution of dimension \(J\). However, if \(J\) is large, this can be computationally infeasible. One advantage of the CAR model structure is that it can be equivalently represented as the joint distribution of \(J\) many univariate normal random variables, where the distribution of each variable is conditional on the remaining \(J-1\) variables. Better still, the distribution of a given \(u_{j,t}\) only depends on it’s adjacent neighbors (hence the name, ‘conditional autoregressive model’), so we can recursively update each \(u_{j,t}\) quite efficiently.

To do this, let \(\mathcal{J}_j\) be the set of neighbor indices of areal unit \(j\). We assume that \(\mathcal{J}_j\) does not change over time. Similarly, let \(d_j\) be the number of neighbors belonging to areal unit \(j\). Then we have

\[\begin{equation}\label{eq:conditional_car} u_{j,t} | u_{k, t, k \neq j} \sim N \left( u_{j, t-1} + \rho_t \sum_{k \in \mathcal{J}_j} \frac{u_{k,t} - u_{k, t-1}}{d_j}, \frac{\tau_t^2}{d_j} \right). \end{equation}\]

The intuition behind equation (\(\ref{eq:conditional_car}\)) is that the distribution of \(u_{j,t}\) given all other \(u_{k, t, k \neq t}\) is normally distributed with a mean that is a weighted average of its neighbors centered on \(u_{j, t-1}\), and with scaled variance which decreases as \(d_j\) increases. Using this form, we can derive the full conditional distribution of each spatio-temporal random effect \(u_{j,t}\) in turn, leading to \(J+T\) many \(1\)-dimensional operations rather than \(T\)-many \(J\)-dimensional operations. In practice, this may be slower for small \(J\), but when \(J\) is large, the computational savings can be substantial.

Writing (\(\ref{eq:posterior}\)) in terms of univariate distributions, we have

\[\begin{equation}\label{eq:uni_posterior} \begin{aligned} &\prod_{t = 1}^T \prod_{j = 1}^J N(y_{j,t} \mid u_{j,t} + \textbf{x}_{j,t} \boldsymbol{\beta}, \sigma^2) \times \prod_{t = 1}^T \prod_{j = 1}^J N \left(u_{j,t} \mid u_{j, t-1} + \rho_t \sum_{k \in \mathcal{J}_j} \frac{u_{k,t} - u_{k, t-1}}{d_j}, \frac{\tau_t^2}{d_j} \right) \\ &\quad \times MVN(\boldsymbol{\beta}\mid \boldsymbol{\mu} _\beta, \boldsymbol{\Sigma} _\beta) \times IG(\sigma^2 \mid a, b) \times \prod_{t = 1}^T IG(\tau^2_t \mid a_t, b_t) \times \prod_{t = 1}^T U(\rho_t \mid c_t, d_t). \end{aligned} \end{equation}\]

We can then see that \(u_{j,t}\) appears in \(3\) terms of (\(\ref{eq:uni_posterior}\)), namely the , the , and the . We do not consider where \(u_{j,t}\) appears in the summation terms of it’s neighbor distributions, because for these distributions the value of \(u_{j,t}\) is being conditioned on, it is treated as fixed.

The product of these three distributions is

Since all three terms are normally distributed, we know that their product will also be normally distributed. Hence, we can find the mean and variance of the resulting full conditional normal distribution by . Thus far, we have discussed incompleting the square for products of multivariate normal distributions, but the same technique will work for products of univariate normals.

Specifically, the resulting full conditional distribution will be distributed as \(N(Vv, V)\), where \(V\) and \(v\) will appear in the product above in the form \(u_{j,t}^2 V^{-1} - 2 u_{j,t} v + Vv^2\). We then have

\[\begin{align*} &\textcolor{Emerald}{N(y_{j,t} \mid u_{j,t} + \textbf{x}_{j,t} \boldsymbol{\beta}, \sigma^2)} \times \textcolor{RoyalBlue}{N \left(u_{j,t} \mid u_{j, t-1} + \rho_t \sum_{k \in \mathcal{J}_j} \frac{u_{k,t} - u_{k, t-1}}{d_j}, \frac{\tau_t^2}{d_j} \right)} \\ &\quad \times \textcolor{OliveGreen}{N \left(u_{j,t+1} \mid u_{j, t} + \rho_{t+1} \sum_{k \in \mathcal{J}_j} \frac{u_{k,t+1} - u_{k, t}}{d_j}, \frac{\tau_{t+1}^2}{d_j} \right)} \\ &\propto \textcolor{Emerald}{\exp \left( \frac{1}{\sigma^2} (y_{j,t} - u_{j,t} - \textbf{x}_{j,t} \boldsymbol{\beta})^2 \right)} \times \textcolor{RoyalBlue}{\exp \left( \frac{d_j}{\tau^2_t} \left( u_{j,t} - u_{j,t-1} - \rho_t \sum_{k \in \mathcal{J}_j} \frac{u_{k,t} - u_{k, t-1}}{d_j} \right)^2 \right)} \\ &\quad \times \textcolor{OliveGreen}{\exp \left( \frac{d_j}{\tau^2_{t+1}} \left( u_{j,t+1} - u_{j,t} - \rho_{t+1} \sum_{k \in \mathcal{J}_j} \frac{u_{k,t+1} - u_{k, t}}{d_j} \right)^2 \right)} \\ &\propto \exp \left( \textcolor{Emerald}{\frac{1}{\sigma^2} (y_{j,t} - u_{j,t} - \textbf{x}_{j,t} \boldsymbol{\beta})^2} + \textcolor{RoyalBlue}{\frac{d_j}{\tau^2_t} \left( u_{j,t} - u_{j,t-1} - \rho_t \sum_{k \in \mathcal{J}_j} \frac{u_{k,t} - u_{k, t-1}}{d_j} \right)^2} \right. \\ &\quad + \left. \textcolor{OliveGreen}{\frac{d_j}{\tau^2_{t+1}} \left( u_{j,t+1} - u_{j,t} - \rho_{t+1} \sum_{k \in \mathcal{J}_j} \frac{u_{k,t+1} - u_{k, t}}{d_j} \right)^2} \right) \\ &\propto \exp \left( \textcolor{Emerald}{\frac{1}{\sigma^2} (u_{j,t}^2 - 2u_{j,t} y_{j,t} + 2 u_{j,t} \textbf{x}_{j,t} \boldsymbol{\beta})} + \textcolor{RoyalBlue}{\frac{d_j}{\tau^2_t} \left( u_{j,t}^2 - 2 u_{j,t} u_{j,t-1} - 2 u_{j,t} \rho_t \sum_{k \in \mathcal{J}_j} \frac{u_{k,t} - u_{k, t-1}}{d_j} \right)} \right. \\ &\quad \left. + \textcolor{OliveGreen}{\frac{d_j}{\tau^2_{t+1}} \left( u_{j,t}^2 - 2 u_{j,t} u_{j,t+1} + 2 u_{j,t} \rho_{t+1} \sum_{k \in \mathcal{J}_j} \frac{u_{k,t+1} - u_{k, t}}{d_j} \right)} \right) \\ &= \exp \left[ u_{j,t}^2 \left( \textcolor{Emerald}{\left( \frac{1}{\sigma^2} \right)} + \textcolor{RoyalBlue}{\left( \frac{d_j}{\tau_t^2} \right)} + \textcolor{OliveGreen}{\left( \frac{d_j}{\tau^2_{t+1}} \right)} \right) \right. \\ &\quad \quad \quad \quad \left. -2 u_{j,t} \left( \textcolor{Emerald}{\left( \frac{1}{\sigma^2} \right) (y_{j,t} - \textbf{x}_{j,t} \boldsymbol{\beta})} + \textcolor{RoyalBlue}{\left( \frac{d_j}{\tau_t^2} \right) \left( u_{j,t-1} + \rho_t \sum_{k \in \mathcal{J}_j} \frac{u_{k,t} - u_{k, t-1}}{d_j} \right)} \right. \right. \\ &\quad \quad \quad \quad + \left. \left. \textcolor{OliveGreen}{\left( \frac{d_j}{\tau_{t+1}^2} \right) \left( u_{j,t+1} - \rho_{t+1} \sum_{k \in \mathcal{J}_j} \frac{u_{k,t+1} - u_{k, t}}{d_j} \right)} \right) + \cdots \right] \end{align*}\]

So we can identify

\[\begin{align*} V^{-1} &= \textcolor{Emerald}{\left( \frac{1}{\sigma^2} \right)} + \textcolor{RoyalBlue}{\left( \frac{d_j}{\tau_t^2} \right)} + \textcolor{OliveGreen}{\left( \frac{d_j}{\tau^2_{t+1}} \right)} \\ v &= \textcolor{Emerald}{\left( \frac{1}{\sigma^2} \right) (y_{j,t} - \textbf{x}_{j,t} \boldsymbol{\beta})} + \textcolor{RoyalBlue}{\left( \frac{d_j}{\tau_t^2} \right) \left( u_{j,t-1} + \rho_t \sum_{k \in \mathcal{J}_j} \frac{u_{k,t} - u_{k, t-1}}{d_j} \right)} \\ &\quad + \textcolor{OliveGreen}{\left( \frac{d_j}{\tau_{t+1}^2} \right) \left( u_{j,t+1} - \rho_{t+1} \sum_{k \in \mathcal{J}_j} \frac{u_{k,t+1} - u_{k, t}}{d_j} \right)} \end{align*}\]

and we can update \(u_{j,t}\) via it’s full conditional as \(N(Vv, v)\).

For \(t = T\), we simply ignore the terms corresponding to the , and we have

\[\begin{align*} V^{-1} &= \textcolor{Emerald}{\left( \frac{1}{\sigma^2} \right)} + \textcolor{RoyalBlue}{\left( \frac{d_j}{\tau_T^2} \right)} \\ v &= \textcolor{Emerald}{\left( \frac{1}{\sigma^2} \right) (y_{j,T} - \textbf{x}_{j,T} \boldsymbol{\beta})} + \textcolor{RoyalBlue}{\left( \frac{d_j}{\tau_T^2} \right) \left( u_{j,T-1} + \rho_T \sum_{k \in \mathcal{J}_j} \frac{u_{k,T} - u_{k, T-1}}{d_j} \right)} \end{align*}\]

and we can update \(u_{j,T}\) via it’s full conditional as \(N(Vv, v)\).

Update \(\sigma^2\)

Looking at (\(\ref{eq:posterior}\)), we see that the term \(\sigma^2\) appears in \(T+1\) terms, namely the and . The product of these terms is takes the form,

\[\begin{equation*} \textcolor{Emerald}{\prod_{t=1}^T MVN(\textbf{y}_t \mid \textbf{u}_t + \textbf{X}_t \boldsymbol{\beta}, \sigma^2 \textbf{I}_J)} \times \textcolor{DeepPink}{IG(\sigma^2 \mid a, b)}. \end{equation*}\]

We know that the inverse gamma distribution (\(\ref{eq:ig}\)) is conjugate for the multivariate normal distribution (\(\ref{eq:mvn}\)), and hence the resulting full conditional distribution will also be distributed as Inverse Gamma. Our goal then is to multiply these densities and identify the resulting expression as an Inverse Gamma distribution for \(\sigma^2\). Specifically,

\[\begin{align*} &\textcolor{Emerald}{\prod_{t=1}^T MVN(\textbf{y}_t \mid \textbf{u}_t + \textbf{X}_t \boldsymbol{\beta}, \sigma^2 \textbf{I}_J)} \times \textcolor{DeepPink}{IG(\sigma^2 \mid a, b)} \\ &= \textcolor{Emerald}{\prod_{t=1}^T (2 \pi)^{-J/2} |\sigma^2 \textbf{I}_J|^{-1/2} \exp\left(-\frac{1}{2} (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta})^\top (\sigma^2 \textbf{I}_J)^{-1} (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta})\right)} \\ &\quad \quad \times \textcolor{DeepPink}{\frac{b^a}{\Gamma (a)} (\sigma^2)^{-a - 1} \exp\left( \frac{-b}{\sigma^2} \right)} \\ &\propto (\sigma^2)^{\textcolor{Emerald}{\frac{-JT}{2}} \textcolor{DeepPink}{-a-1}} \exp \left( \frac{\textcolor{Emerald}{-\frac{1}{2} \sum_{t=1}^T (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta})^\top (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta})} \textcolor{DeepPink}{-b}}{\sigma^2}\right) \end{align*}\]

which we can identify as in the form of an Inverse Gamma pdf (\(\ref{eq:ig}\)), with parameters

\[\begin{align*} \tilde{a} &= \textcolor{DeepPink}{a} + \textcolor{Emerald}{\frac{JT}{2}} \\ \tilde{b} &= \textcolor{DeepPink}{b} + \textcolor{Emerald}{\frac{1}{2} \sum_{t=1}^T (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta})^\top (\textbf{y}_t - \textbf{u}_t - \textbf{X}_t \boldsymbol{\beta})} \end{align*}\]

and we can update \(\sigma^2\) via it’s full conditional as \(IG(\tilde{a}, \tilde{b})\).

Update \(\tau^2_t, \quad t = 1, \ldots, T\)

When viewing the full posterior (\(\ref{eq:posterior}\)), we see that \(\tau_t^2\) appears in two terms, namely in and in , which is distributed as inverse gamma. The product of these distributions will be

\[\begin{equation*} \textcolor{RoyalBlue}{MVN(\textbf{u}_t \mid \textbf{u}_{t-1}, \tau_t^2 \textbf{Q}(\rho_t))} \times \textcolor{DarkMagenta}{IG(\tau_t^2 \mid a_t, b_t)} \end{equation*}\]

which will form a Inverse Gamma full conditional distribution by conjugacy. As before, we multiply these densities to obtain

\[\begin{align*} &\textcolor{RoyalBlue}{MVN(\textbf{u}_t \mid \textbf{u}_{t-1}, \tau_t^2 \textbf{Q}(\rho_t))} \times \textcolor{DarkMagenta}{IG(\tau_t^2 \mid a_t, b_t)} \\ &= \textcolor{RoyalBlue}{(2 \pi)^{-J/2} |\tau_t^2 \textbf{Q}(\rho_t)|^{-1/2} \exp \left( -\frac{1}{2} (\textbf{u}_t - \textbf{u}_{t-1})^\top (\tau_t^2 \textbf{Q}(\rho_t))^{-1} (\textbf{u}_t - \textbf{u}_{t-1}) \right)} \\ &\quad \quad \times \textcolor{DarkMagenta}{\frac{b_t^{a_t}}{\Gamma (a_t)} (\tau_t^2)^{-a_t - 1} \exp \left( \frac{-b_t}{\tau_t^2} \right)} \\ &\propto (\tau_t^2)^{\textcolor{RoyalBlue}{- \frac{J}{2}} \textcolor{DarkMagenta}{-a_t-1}} \exp \left( \frac{\textcolor{RoyalBlue}{- \frac{1}{2} (\textbf{u}_t - \textbf{u}_{t-1})^\top (\textbf{u}_t - \textbf{u}_{t-1})} \textcolor{DarkMagenta}{-b_t}}{\tau_t^2}\right) \end{align*}\]

which we can identify as in the form of an Inverse Gamma pdf (\(\ref{eq:ig}\)), with parameters

\[\begin{align*} \tilde{a_t} &= \textcolor{DarkMagenta}{a_t} + \textcolor{RoyalBlue}{\frac{J}{2}} \\ \tilde{b_t} &= \textcolor{DarkMagenta}{b_t} + \textcolor{RoyalBlue}{\frac{1}{2} (\textbf{u}_t - \textbf{u}_{t-1})^\top (\textbf{u}_t - \textbf{u}_{t-1})} \end{align*}\]

and we can update \(\tau_t^2\) via it’s full conditional as \(IG(\tilde{a_t}, \tilde{b_t})\).

Update \(\rho_t, \quad t = 1, \ldots, T\)

Until now we have worked directly with full conditional distributions for model parameters via conjugate priors. This has made things easier, as we could directly determine the form of the full conditional distribution. Now, we introduce a to update \(\rho_t, t = 1, \ldots, T\). The metropolis step is essentially a process whereby you propose a new value for your parameter, and accept or reject that proposed value based on a proposal distribution. In (\(\ref{eq:mod}\)), we chose a Uniform prior for \(\rho_t\) as \(U(c_t, d_t)\). Importantly, the Uniform distribution has bounded support, meaning our proposal distribution (which is typically Guassian) may propose values outside the support of \(\rho_t\). Hence, we will introduce a for \(\rho_t\), which will give it support on the real line \(\mathbb{R}\). This way, we can use a Gaussian proposal distribution, evaluate the acceptance probability of this transformed variable value, and then back transform our sample upon acceptance. We will outline the general framework here.

First, consider the full posterior (\(\ref{eq:posterior}\)). As we have done in all steps, we identify the terms in which \(\rho_t\) appears, namely in and in , which is Uniform on \((c_t, d_t)\). The product of these densities is

\[\begin{align*} &\textcolor{RoyalBlue}{MVN(\textbf{u}_t \mid \textbf{u}_{t-1}, \tau_t^2 \textbf{Q}(\rho_t))} \times \textcolor{Chocolate}{U(\rho_t \mid c_t, d_t)} \\ &= \textcolor{RoyalBlue}{(2 \pi)^{-J/2} |\tau_t^2 \textbf{Q}(\rho_t)|^{-1/2} \exp \left( - \frac{1}{2} (\textbf{u}_t - \textbf{u}_{t-1})^\top (\tau_t^2 \textbf{Q}(\rho_t))^{-1} (\textbf{u}_t - \textbf{u}_{t-1}) \right)} \\ &\quad \quad \times \textcolor{Chocolate}{\frac{1}{d_t - c_t} \boldsymbol{1}_{c_t < \rho_t < d_t}} \end{align*}\]

where \(\boldsymbol{1}\) is the indicator function. Next, we employ a transformation on the Uniform variable \(\rho_t\) of the form

\[\begin{equation*} g(\rho_t) = \log \left( \frac{\rho_t - c_t}{d_t - \rho_t} \right). \end{equation*}\]

To get the distribution of this newly transformed variable \(g(\rho_t)\), we require the Jacobian. First, we find the inverse function of the transformation, which is

\[\begin{equation*} g^{-1}(g(\rho_t)) = \frac{d_t e^{g(\rho_t)} + c_t}{e^{g(\rho_t)} + 1}. \end{equation*}\]

Next, we take the derivative of this function with respect to \(g(\rho_t)\), which gives

\[\begin{equation*} \frac{e^{g(\rho_t)} (d_t - c_t)}{(e^{g(\rho_t)} + 1)^2}, \end{equation*}\]

and plugging the expression for \(g(\rho_t)\) into this equation yields

\[\begin{equation*} \frac{\frac{\rho_t - c_t}{d_t - \rho_t} (d_t - c_t)}{\left(\frac{\rho_t - c_t}{d_t - \rho_t} + 1 \right)^2} \propto \frac{\frac{\rho_t - c_t}{d_t - \rho_t}}{\left(\frac{d_t - \rho_t + \rho_t - c_t}{d_t - \rho_t} \right)^2} \propto \frac{\frac{\rho_t - c_t}{d_t - \rho_t}}{\frac{1}{(d_t - \rho_t)^2}} = (\rho_t - c_t)(d_t - \rho_t). \end{equation*}\]

With this jacobian adjustment, we can now determine the to which we will add the Jacobian adjustment. We have

\[\begin{align*} &\log \left[ \textcolor{RoyalBlue}{(2 \pi)^{-J/2} |\tau_t^2 \textbf{Q}(\rho_t)|^{-1/2} \exp \left( - \frac{1}{2} (\textbf{u}_t - \textbf{u}_{t-1})^\top (\tau_t^2 \textbf{Q}(\rho_t))^{-1} (\textbf{u}_t - \textbf{u}_{t-1}) \right)} \right. \\ &\quad \quad \left. \times \textcolor{Chocolate}{\frac{1}{d_t - c_t} \boldsymbol{1}_{c_t < \rho_t < d_t}} \times (\rho_t - c_t)(d_t - \rho_t) \right] \\ &\propto \textcolor{RoyalBlue}{- \frac{1}{2} \log |\tau_t^2 \textbf{Q}(\rho_t)| - \frac{1}{2} (\textbf{u}_t - \textbf{u}_{t-1})^\top (\tau_t^2 \textbf{Q}(\rho_t))^{-1} (\textbf{u}_t - \textbf{u}_{t-1})} + \log (\rho_t - c_t) + \log (d_t - \rho_t) \end{align*}\]

We can then proceed by generating a proposed value of \(g(\rho_t)\) from some proposal distribution (usually Gaussian), evaluating the log target density plus Jacobian adjustment, and then accepting or rejecting the back-transformed \(g^{-1}(g(\rho_t))\) based on this acceptance probability.