Hi all,

I am a complete beginner to Stan, so apologies in advance if I have made an obvious mistake or asking a silly question.

**Set up**

I am attempting to fit a “rolling” regression model in which the slope parameter varies through time:

y_t \sim \text{normal}(\mu_t, \sigma)

\mu_t = \alpha + \beta_t\cdot X_t

In particular, I want to apply a Gaussian process prior to \beta so that it varies relatively smoothly through time:

\beta \sim \text{multi-normal}\left(0, k(t_i, t_j) \right)

Where k is an exponentiated quadratic kernel.

**Problem**

When I use a fixed (true) GP as the prior, and infer \beta_t, \alpha, and \sigma, I can get back a posterior that closely captures the true values. However, when I try to also infer one of the GP parameters, even with an informative prior, sampling never seems to progress past the first iteration / warmup. Specifically I get stuck at: “Iteration: 1 / 1000 [ 0%] (Warmup)”.

Ultimately of course I would like to infer all parameters and build up a more complicated model. I understand that modelling \beta as a random walk is another option. I have done this already and would like to trouble-shoot the GP prior for now.

```
# Simulation parameters
gp_alpha <- 0.25
gp_rho <- 20
alpha <- 0
sigma <- 0.05
# Set up GP
eps <- sqrt(.Machine$double.eps)
t <- 1:500
D <- as.matrix(dist(t))^2
K <- gp_alpha^2 * exp(-1/(2 * gp_rho^2) * D) + diag(eps, length(t))
# Simulate data
set.seed(9)
beta <- mvtnorm::rmvnorm(1, mean=rep(0, length(t)), sigma=K)[1, ] # unobserved
x <- runif(length(t)) # observed
y <- alpha + (beta * x) + rnorm(length(t), sd=sigma) # observed
# plot(t, beta, type = 'l', lwd=2, xlab="time", ylab="Beta")
```

```
data {
int<lower=1> N;
vector[N] y;
real x[N];
real t[N];
real<lower=0> gp_rho;
}
parameters {
real alpha;
real<lower=0> sigma;
real<lower=0> gp_alpha;
vector[N] Z;
}
model {
// Partially-fixed GP prior
matrix[N, N] cov = cov_exp_quad(t, gp_alpha, gp_rho)
+ diag_matrix(rep_vector(1e-10, N));
matrix[N, N] L_cov = cholesky_decompose(cov);
vector[N] beta = L_cov * Z;
// Mean model
for (i in 1:N) {
y[i] ~ normal(alpha + beta[i] * x[i], sigma);
}
Z ~ normal(0, 1);
// other priors
gp_alpha ~ normal(0, 0.2);
sigma ~ exponential(1);
alpha ~ normal(0, 0.05);
}
```

```
# Fit model
dat <- list(N = length(y), y = y, x = x, t = 1:length(y), gp_rho = gp_rho)
s <- rstan::stan_model(model_code = mod, verbose = TRUE)
fit <- rstan::sampling(s, data = dat, cores = 4, iter = 1000, chains = 4)
```

Any help or advice is appreciated. I am using R version 4.1.0 with rstan version 2.21.2.