I am sharing the description of the model along with the data and code to run the model. I am getting error while fitting the model. Kindly suggest possible solution for this problem. Thanks.
Note: I have uploaded the csv file of data here. However in code, it is written for xlsx file. Please modify the code to import csv file of inflation data.
Model description.pdf (58.5 KB)
01_Write Model Fertility.R (7.1 KB)
02_Data Prep.R (3.7 KB)
03_HB Forecast.R (3.9 KB)
Inflation-data.csv (93.9 KB)
Hi, @Rahul and welcome to the Stan forums. And thanks for including a reproducible example.
Your model description is relatively straightforward (the missing close paren in the first normal threw me a bit). I’m not exactly sure why it took so much Stan code, but maybe there are just a lot of different variables floating around.
I would strongly recommend writing a standalone Stan model rather than embedding it in R code. With a standalone version, you can get line numbers for errors in a meaningful way. Here’s the Stan code from your R script:
data {
int<lower=0> A; // number of age groups
int<lower=0> Y; // number of time periods
int<lower=0> L; // number of locations
int<lower=0> R; // number of regions
int<lower=0> S; // number of super-regions
int<lower=1,upper=R> region[L]; // region for each location
int<lower=1,upper=S> super_region[L]; // super-region for each location
int<lower=1,upper=S> super_region_r[R]; // super-region for each region
real<lower=0,upper=1> f[A, Y, L]; // fertility rate
real z[A, Y, L]; // predictor variable Z for each location and time
real X[A, Y, L]; // age data for each location, age group, and time
}
transformed data {
// Create a 3D array for transformed x
real Z[A, Y, L];
// Flatten the 3D array z into a 1D array to apply min
real z_flat[A * Y * L];
int index = 1;
for (l in 1:L) {
for (y in 1:Y) {
for (a in 1:A) {
z_flat[index] = z[a, y, l];
index = index + 1;
}
}
}
// Find the minimum value across all elements of z_flat
real shift_value = abs(min(z_flat)); // Find the most negative value in z and shift to zero
// Adjust the x values by adding the shift value
for (l in 1:L) {
for (y in 1:Y) {
for (a in 1:A) {
Z[a, y, l] = z[a, y, l] + shift_value;
}
}
}
}
parameters {
// Random intercepts and slopes
real alpha[A, Y, L]; // Random intercepts for each location, age, and time
real kappa[A, Y, L]; // Random slopes for each location, age, and time
real beta[A, Y, L]; // Random slopes for Z with a two-year lag
real alpha_mu; // Mean for the random intercept
real<lower=0> sigma_f; // Standard deviation for the fertility model
real<lower=0> sigma_alpha;
real<lower=0> sigma_kappa; // Standard deviation for the second-order random walk for kappa
real<lower=0> sigma_beta;
// Temporal effects
real theta[Y, L]; // Temporal effects for each location and time
real phi[L]; // AR(1) coefficients for each location
real<lower=0> sigma_theta[L]; // Standard deviations for temporal error term for each location
// Location effects across different levels of hierarchy
real eta[L]; // location-specific effects
real eta_r[R]; // region-specific effects
real eta_sr[S]; // super-region-specific effects
real eta_w; // world-wide effect
real<lower=0> sigma_eta_l; // Standard deviation of location effects
real<lower=0> sigma_eta_r; // Standard deviation of regional effects
real<lower=0> sigma_eta_sr; // Standard deviation of super-region effects
// Residuals for ARIMA model (epsilon_a,y,l)
real epsilon[A, Y, L]; // residuals for location-age-time
real phi_epsilon[S]; // AR(1) coefficients for each super-region
real<lower=0> sigma_epsilon[S]; // Standard deviations for residuals in super-region s
}
model {
sigma_f ~ normal(0, 1);
alpha_mu ~ normal(0, 10);
// Priors for AR(1) coefficient and residual standard deviation for each super-region
sigma_alpha ~ normal(0, 1);
sigma_beta ~ normal(0, 1);
sigma_kappa ~ normal(0, 1); // Prior for standard deviation of the random walk
phi ~ normal(0, 0.5);
sigma_theta ~ normal(0, 1);
sigma_eta_l ~ normal(0, 1);
sigma_eta_r ~ normal(0, 1);
sigma_eta_sr ~ normal(0, 1);
phi_epsilon ~ normal(0, 0.5);
sigma_epsilon ~ normal(0, 1);
// Random intercept alpha
for (l in 1:L) {
for (y in 1:Y) {
for (a in 1:A) {
alpha[a, y, l] ~ normal(alpha_mu, sigma_alpha);
}
}
}
// Second-order random walk (RW2) for kappa over age
for (l in 1:L) {
for (y in 1:Y) {
// Initial condition for first age group
kappa[1, y, l] ~ normal(0, sigma_kappa);
kappa[2, y, l] ~ normal(kappa[1, y, l], sigma_kappa);
for (a in 3:A) {
// RW2 process: second-order random walk for kappa
kappa[a, y, l] ~ normal(2 * kappa[a-1, y, l] - kappa[a-2, y, l], sigma_kappa);
}
}
}
// Beta for lagged Z (t-2) effect
for (l in 1:L) {
for (y in 1:Y) {
for (a in 1:A) {
beta[a, y, l] ~ normal(0, sigma_beta); // Priors for beta (location, age, time-specific)
}
}
}
// Temporal effects (AR(1) structure)
for (l in 1:L) {
theta[1,l] ~ normal(0, sigma_theta[l]); // Initial condition for time 1
for (y in 2:Y) {
theta[y, l] ~ normal(phi[l] * theta[y-1, l], sigma_theta[l]); // AR(1) process
}
}
// location-specific effects (eta_l)
// Level 1: eta_l ~ N(eta_r[l], sigma_eta_l)
for (l in 1:L) {
eta[l] ~ normal(eta_r[region[l]], sigma_eta_l); // eta_l depends on eta_r[region[l]]
}
// Level 2: eta_r ~ N(eta_sr[r], sigma_eta_r)
for (r in 1:R) {
eta_r[r] ~ normal(eta_sr[super_region_r[r]], sigma_eta_r); // eta_r depends on eta_sr[super_region[r]]
}
// Level 3: eta_sr ~ N(eta_w, sigma_eta_sr)
for (s in 1:S) {
eta_sr[s] ~ normal(eta_w, sigma_eta_sr); // eta_sr depends on eta_w
}
// Level 4: eta_w ~ N(0, 1)
eta_w ~ normal(0, 1); // eta_w is drawn from N(0, 1)
// Residuals (epsilon_l,a,t) modeled with ARIMA(1,0,0) for each super-region
for (s in 1:S) { // Loop over super-regions
// Loop through locations within super-region s
for (l in 1:L) {
if (super_region[l] == s) { // Locations in super-region s
for (a in 1:A) {
// Prior for epsilon at t = 1
epsilon[a, 1, l] ~ normal(0, sigma_epsilon[super_region[l]]); // Residual for time 1
for (y in 2:Y) {
epsilon[a, y, l] ~ normal(phi_epsilon[super_region[l]] * epsilon[a, y-1, l], sigma_epsilon[super_region[l]]);
}
}
}
}
}
// Likelihood
real logit_f_star[A,Y,L];
for (l in 1:L) {
for (y in 1:Y) {
for (a in 1:A) {
real logit_f = alpha[a, y, l] + kappa[a, y, l] * X[a, y, l]
+ beta[a, y, l] * Z[a, y, l] + theta[y, l] + eta[l] + epsilon[a, y, l];
logit_f_star[a, y, l] ~ normal(logit_f, sigma_f);
// Logistic transformation for observed fertility
f[a, y, l] ~ normal(inv_logit(logit_f_star[a, y, l]), sigma_f);
}
}
}
}
generated quantities {
real f_pred[A, Y, L]; // Predicted fertility rates (as a 3D array)
// Generate predicted values
for (l in 1:L) {
for (y in 1:Y) {
for (a in 1:A) {
// Generate predicted values for fertility with logit transformation
real logit_f = alpha[a, y, l] + kappa[a, y, l] * X[a, y, l]
+ beta[a, y, l] * Z[a, y, l] + theta[y, l] + eta[l]
+ epsilon[a, y, l];
f_pred[a, y, l] = inv_logit(normal_rng(logit_f, sigma_f));
}
}
}
}
This is a huge model that’s way too much for someone to debug on the forums. What I would recommend is trying to build this model up from scratch one piece at a time. What is the most complicated model you were able to run before things fell apart? If the answer is none, I would suggest trying to strip this model down to its components and then build back up. then you can see where things go wrong, whereas it’s almost impossible to debug a Stan model of this size if you don’t have the next simpler model already working.
there are things in your model that look very unusual and I suspect are wrong, like this:
f[a, y, l] ~ normal(inv_logit(logit_f_star[a, y, l]), ...
It’s very unusual to restrict the location of an unconstrained distribution like the normal to be in the range (0, 1). I see that in your data generation in generated quantities, the logit’s on the outside:
f_pred[a, y, l] = inv_logit(normal_rng(logit_f, sigma_f));