Modelling the heteroscedasticity (similar to varExp)

Modeling
1

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 thatimage

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

2

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