Modelling the heteroscedasticity (similar to varExp)

Hi,

I would like become familiar with Bayesian methods but I don’t have many opportunities. Now, I have worked with a model (frequentist approach) and I am trying to replicate it in rstan. Because my data suffer from heteroscedasticity I have modeled that too by taking the power of the covariant which lookls like that

As example let’s consider the following data

``````library(rstan)
library(nlme)

set.seed(1234)
x1 <- rep(1:100, 2)
x2 <- runif(200, 10, 100)
summary(x2)

a <- 0
b1 <- 1
b2 <- 1
sigma2 <- x1^1.3
eps <- rnorm(200, mean = 0, sd = sqrt(sigma2))
y <- a + b1*x1 + b2*x2 + eps
dat <- data.frame(x1, x2, y)
``````

The model I run look like this

``````VMat <- varPower(form = ~x1)
mod <- gls(y ~ x1 + x2, weights = VMat, data = dat)
summary(mod)
``````

Now, I have hard times to replicate that on Stan. For time being the following run smoothly without modeling the heteroscedasticity.

``````

model_stan_lm <- '
data{
int <lower = 1> N;
vector[N] x1;
vector[N] x1_abs;
vector[N] x2;
vector[N] y;

}

parameters {
real a;
real b1;
real b2;
real delta;
real <lower = 0> sigma;
}

transformed parameters {
vector[N] mu;
mu = a + b1*x1 + b2*x2;

//vector[N] var_exp;
//vector[N] sd_exp;
//var_exp = sqrt(sigma) * pow(abs_x1, 2*delta); //covariate variance
//sd_exp = pow(var_exp, 2);
}

model{
y ~ normal(mu, sigma);
}
'

stan_samples_lm <- stan(model_code = model_stan_lm , model_name='lm_trial', verbose = FALSE,
data = list(
N = nrow(dat),
x1 = dat\$x1,
x2 = dat\$x2,
y = dat\$y,
x1_abs = abs(dat\$x1)
)
)
stan_samples_lm
``````

However, when I try to model the heteroscedasticity I am getting problems with the types of data I use but by trying to resolve them I end up in new problems.

``````

model_stan_lm <- '
data{
int <lower = 1> N;
vector[N] x1;
vector[N] x1_abs;
vector[N] x2;
vector[N] y;

}

parameters {
real a;
real b1;
real b2;
real delta;
real <lower = 0> sigma;
}

transformed parameters {
vector[N] mu;

vector[N] var_exp;
vector[N] sd_exp;
mu = a + b1*x1 + b2*x2;

var_exp = sqrt(sigma) * pow(x1_abs, 2*delta);
sd_exp = pow(var_exp, 2);
}

model{
y ~ normal(mu, sd_exp);
}
'

stan_samples_lm <- stan(model_code = model_stan_lm , model_name='lm_trial', verbose = FALSE,
data = list(
N = nrow(dat),
x1 = dat\$x1,
x2 = dat\$x2,
y = dat\$y,
x1_abs = abs(dat\$x1)
)
)
stan_samples_lm
``````

Any Idea how to model the heteroscedasticity;

Thanks

Is it fair to write your generating process:

``````y = rnorm(200, a + b1 * x1 + b2 * x2, sqrt(x1^d));
``````

And you want to estimate a, b1, b2, and d?

If so, then this should be how to code that in Stan:

``````y ~ normal(a + b1 * x1 + b2 * x2, sqrt(x1^d));
``````
1 Like