model <- 'data { int N; int K; int T; matrix[N,K] X[T] ; vector[N * T] y; matrix[N,N] I; } parameters { vector[K] beta; real sigma_nug; real phi; } transformed parameters { vector[N] Y[T]; vector[N] epsilon[T]; // error term vector[N] mu[T]; // mean real var_nug; // nugget for (t in 1:T){; // converts y to a 3D array Y[t] = y[((t - 1) * N + 1):(t * N)]; } var_nug = sigma_nug ^ 2; mu[1] = X[1] * beta; epsilon[1] = Y[1] - mu[1]; for (t in 2:T){ mu[t] = X[t] * beta; epsilon[t] = Y[t] - mu[t]; mu[t] = mu[t] + phi * epsilon[t-1]; } } model { for (t in 1:T){ target += multi_normal_cholesky_lpdf(Y[t] | mu[t], cholesky_decompose(var_nug * I) ); } sigma_nug ~ uniform(0,50); phi ~ uniform(-1, 1); // } generated quantities { vector[N] log_lik[T]; for (t in 1:T) { for (n in 1:N) { log_lik[t][n] = multi_normal_cholesky_lpdf(Y[t] | mu[t], cholesky_decompose(var_nug * I) ) ; } } } '