we are two Italian statistical student of Sapienza Università di Roma. We are working about Testing Hypotheses via a mixture estimation model, based on the work of Kamary, Mergensen, Robert and Rousseau.
In particular we would like to implement this procedure for a mixture of two autoregressive processes AR(1) and AR(2). Our aim is to estimate the weights of the mixture and analyze if their behaviour is coherent with the posterior probability of the model computed through the Bayes Factor.
We tried to implement ourself an Metropolis within Gibbs algorithm for all the parameters of the mixture and it seems to work quite well when the true model is AR(1) but not when the true model is AR(2).
Furthermore we would like to implement the same procedure with the package Rstan. But we had a lot of problems in the code.
Could you help us with same advises on the implementation of our code.
We joint our code in pdf (r-markdown) with a short summary, an attempt of the implementation with Rstan package and the paper based on.
library(rstan)
y_2 = arima.sim(list(order=c(1,0,0),ar = c(0.5)) ,n=10000)
mix <- '
data {
int<lower=0> N;
real y[N];
}
parameters {
real<lower=0,upper=1> theta;
real alpha;
real mu;
real beta;
real gamma;
real<lower=0> sigma;
real<lower=0> tau;
}
model {
theta ~ beta(1,1);
mu ~ normal(0.5,6);
sigma ~ inv_gamma(2, 0.1);
for(n in 3:N)
target += log_mix(theta, normal_lpdf(y[n] | mu + beta * y[n-1], sigma),
normal_lpdf(y[n] |alpha + beta * y[n-1] + gamma * y[n-2], tau));
}
'
dat = list(y_2)
test = stan(model_code = mix, model_name = "example",
data = dat, iter = 10000, chains = 2,
verbose = F)
(Just FYI, I added a code fence around the code in the original post. To do that, it’s ``` to start and ``` to end)
Taking a quick look at the Stan code, it looks like either there’s an error in the model or I’m not understanding something about your model. In the log_mix() line, you’ve included the same beta in both terms. That doesn’t seem right.
I’d suggest starting a little simpler: write the model with just the AR(1) term and then write the AR(2) model to make sure they both work. It looks like you don’t have priors on all the parameters; it’ll depend on the data whether or not the parameters are identified.
[this is duplicated from my response to the original]
I’m not familiar with their comparison approach, but I keep seeing it on Robert’s blog, so I should probably get up to speed.
Are you sure you want to share that beta parameter across both models being compared? We don’t recommend inverse gamma priors and you only have a prior on sigma, not tau. You also don’t have priors on man of the other parameters.
Did you choose a set of parameters for simulation that ensures stationarity of the series?
Is the ARIMA class of models the same as the AR models?
P.S. Indentation really helps in reading programs, but maybe it just got removed when it wasn’t escaped for code. I added the escapes.
P.P.S. That’s “Mengersen” and “Kamary”. At least it’s open access (go arXiv).