data {
  int<lower=0> n;
  real y[n];
  real<lower=0> sigma1;
}

parameters {
  real mu;
  real<lower=0> alpha0;
  real<lower=0,upper=1> alpha1;
  real<lower=0,upper=(1-alpha1)> beta1;
}

transformed parameters {
  real<lower=0> sigma[n];
  sigma[1] = sigma1;
  for (t in 2:n)
    sigma[t] = sqrt(alpha0
                     + alpha1 * pow(y[t-1] - mu, 2)
                     + beta1 * pow(sigma[t-1], 2));
}

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