I’m trying to port a hierarchical repeat sales index model from BUGS , the original code from the paper is attached (I think the author is @avdminne, ) . My Stan model compiles and fits, but I need to turn max_treedepth up to at least 14 in order to avoid warnings - usually all or almost all iterations exceed the treedepth at default settings. This is fine but taking max_treedepth from 10 to 14 lengthens the runtime 6x+ and I’m already only using a subset of the data. Is there anything I could try to avoid a higher treedepth, maybe reparameterization tricks or alternative priors?
I also read on the “Guide to Errors” about looking for parameters correlated with energy_ , but I have a lot of parameters and since all iterations typically exceed the max_treedepth I wonder whether I’m defining something in a not ideal way.
data {
int<lower=1> N; // N observations
int<lower=1> Na; // N areas
int<lower=1> Nt; //N timeperiods
vector[N] y; // log price diff
int<lower=1> month1[N]; //month of repeat trxn
int<lower=1> month0[N]; //month of first trxn
int<lower=1> area_id[N]; // area id
real<lower=0> adf; // parameter for DoF
}
parameters {
real<lower=0> sigma_y;
real beta[Nt];
real inc[Nt];
real kappa[Nt];
real<lower=0> tau_alpha;
real<lower=0> tau_kappa;
real<lower=0> tau_eta;
matrix[Na,Nt] alpha;
real<lower=1> nu; // degrees of freedom for t
}
transformed parameters {
real mu_beta[Nt];
real mu[N];
mu_beta[1] = 0;
for (t in 2:Nt){
mu_beta[t] = beta[t-1]+inc[t]*kappa[t]; // 'matt trick?'
}
for (i in 1:N){
mu[i] = (beta[month1[i]] - beta[month0[i]]) +
(alpha[area_id[i], month1[i]] - alpha[area_id[i], month0[i]]);
}
}
model {
// Priors
sigma_y ~ inv_gamma(0.1,0.1);
tau_alpha ~ lognormal(5,1);
tau_kappa ~ inv_gamma(0.001, 0.001);
tau_eta ~ inv_gamma(0.001, 0.001);
nu ~ exponential(adf);
beta ~ normal(mu_beta, tau_eta);
inc ~ normal(0,1);
for (t in 2:Nt) {
kappa[t] ~ normal(kappa[t-1], tau_kappa);
}
for (a in 1:Na){
for (t in 2:Nt){
alpha[a,t] ~ normal(alpha[a,t-1], tau_alpha);
}
}
y ~ student_t(nu, mu, sigma_y);
}
I’m omitting the other hierarchical features for property types that is in the BUGS model. Also I tried to incorporate the “Matt trick” as described in this other post Structural Time Series research using Stan.
EDIT: Based on reading discussions elsewhere and the guide to priors it seems like inverse gamma priors are necessary for BUGS but not ideal for Stan, I’ve changed them to standard student_t with DoF between 2-4.
BUGScode-2017-hiearchichal-repeat-sales.pdf (56.1 KB)
This suggests there is some correlation in your parameters. Looking at the model, this is where I would start. try, for a few, i plotting pairs of:
- beta[i], inc[i], and kappa[i]
- mu_beta, tau_eta, and beta[i]
- kappa[i] and kappa[i-1]
- kappa[i] and tau_kappa
- nu and sigma
Also, I’m not totally sure I follow what you are trying to do, but this paper may be of interest, which models repeat sales in Stan.
1 Like
Thanks, the model name is sort-of confusing and I could’ve included some more info. Repeat Sales models in real estate try to estimate price appreciation/depreciation using the log price change between sales for properties that have changed hands multiple times to try and and avoid having to control for hedonic factors of each property.
There are lots of problems with this model, starting with priors that we don’t recommend on either statistical or computational grounds (we have a Wiki on prior recommendations and also the regression chapter in the manual).
The problem may be with the parameterization of beta, which I don’t quite follow. There’s this in the model block
beta ~ normal(mu_beta, tau_eta);
where mu_beta is a transformed parameter and tau_eta a parameter. You want to use the non-centered paramterization on beta, and make thatbeta_raw ~ normal(0, 1);and then definebeta = mu_beta + beta_raw * tau_eta;. That naming scheme is confusing. We'd usually usesigma_betaas the name of the parameter rather thantau_eta`.
Also, do you realize stan uses standard deviation, not variance?
Then you may also need a non-centered parameerization of the time series, but I’d start with beta and by tighteing the priors.
I blindly copied the priors from the original BUGS model without realizing they are on inverse-variance, Ive changed them all based on the wiki guidelines. Sorry for the confusing tau names, I did it to mostly preserve the syntax of the BUGS model and make comparison a little easier. I can’t post the BUGS model in text because it’s trapped in a pdf, but I attached it to the original post.
Anyway, I did the non-centered parameterization and modified all the priors, see below. I’m still getting lots of (but less than before) iterations that exceed max treedepth. Considering this model presumably worked well in BUGS could there be something that makes it unsuitable for use in Stan?
Thanks,
data {
int<lower=1> N; // N observations
int<lower=1> Na; // N areas
int<lower=1> Nt; //N timeperiods
vector[N] y; // log price diff
int<lower=1> month1[N]; //month of repeat trxn
int<lower=1> month0[N]; //month of first trxn
int<lower=1> area_id[N]; // area id
real<lower=0> adf; // parameter for DoF
}
parameters {
real<lower=0> sigma_y;
real beta_raw[Nt];
real inc[Nt];
real kappa[Nt];
real<lower=0> tau_alpha;
real<lower=0> tau_kappa;
real<lower=0> tau_eta;
matrix[Na,Nt] alpha;
real<lower=1> nu; // degrees of freedom for t
}
transformed parameters {
real beta[Nt];
real mu_beta[Nt];
real mu[N];
beta[1] = 0;
mu_beta[1] = 0;
for (t in 2:Nt){
mu_beta[t] = beta[t-1]+inc[t]*kappa[t]; // 'matt trick?'
beta[t] = mu_beta[t] + beta_raw[t] * tau_eta; // non-centered reparameterization
}
for (i in 1:N){
mu[i] = (beta[month1[i]] - beta[month0[i]]) +
(alpha[area_id[i], month1[i]] - alpha[area_id[i], month0[i]]);
}
}
model {
// Priors
sigma_y ~ student_t(2,0,10);
tau_alpha ~ normal(0,10);
tau_kappa ~ normal(0,10);
tau_eta ~ normal(0,10);
nu ~ exponential(adf);
beta_raw ~ normal(0,1);
inc ~ normal(0,1);
for (t in 2:Nt) {
kappa[t] ~ normal(kappa[t-1], tau_kappa);
}
for (a in 1:Na){
for (t in 2:Nt){
alpha[a,t] ~ normal(alpha[a,t-1], tau_alpha);
}
}
y ~ student_t(nu, mu, sigma_y);
}
Could you post the pairs plot for sigma_y, tau_alpha, tau_kappa, tau_eta, nu, beta_raw[i], inc[i], kappa[i], alpha[j,i] for a a few i,j?
The model may not actually have worked well in BUGS. BUGS has fewer diagnostics than Stan, so many models that appear to work in BUGS and failed in Stan, actually were never working in BUGS to begin with.
i didn’t see a link.
can you cut and paste from the pdf?
kappa still has a centered parameterization. do you get divergences? does it work with simulated data?
I basically never get divergences, just exceeding treedepth. Haven’t tried simulated data yet.
Here’s the BUGS code. Note I’m not including delta in Stan.
Here is a pairs plot, I’m not sure what exactly to look for but it seems like alpha being bimodal is worrisome? And alpha is correlated with energy. Pairs for other (i,j) look similar, although sometimes alpha will be unimodal or sometimes even trimodal.
Thanks,
What is alpha in the pairs plot? In your model you don’t have a scalar quantity alpha.
Alpha is an [i,j] matrix, I am plotting just one representative i and j here.
Thanks.
Try using the ‘Matt trick’ on the alphas. Parameterize the time series in terms of unit normal innovation terms, and then get the alphas in the transformed parameter block.
That is,
parameters {
matrix[Na, Nt] alpha_std;
transformed parameters {
matrix[Na, Nt] alpha;
for (a in 1:Na) {
// implicit: alpha[a, 1] ~ normal(mu_alpha0, tau_alpha0)
alpha[a, 1] = mu_alpha0 + tau_alpha0 * alpha_std[a, 1];
for (t in 2:Nt)
// implicit: alpha[a, t] ~ normal(alpha[a, t - 1], tau_alpha)
alpha[a, t] = alpha[a, t - 1] + tau_alpha * alpha_std[a, t - 1]
}
model {
to_vector(alpha_std) ~ normal(0, 1);
It would be a bit more efficient to compute the alpha term by scaling all
the alpha_std at once with a scalar-matrix multiply and then putting alpha together using cumulative_sum.
3 Likes
Thanks, that seems to have solved the problem, no more iterations exceeding max_treedepth as long as I set the max to 11.
I am having a similar problem in converting BUGS code into STAN. Could you upload a full version of the hierarchical trend repeat sales model? I could follow the piecemeal inputs. It will be helpful a lot.
Tagging @sar2160 to make sure they are the post (and hopefully can share the Stan code)
I am trying to solve the same issue but not very successful. Could upload the revised and working Stan code for the Hierarchical Trend Repeat Sales Model with Local Linear Trend and t-distribution? It is very difficult to follow the piecemeal codes.
