Fit a univariate GEV distribution in stan

gen_extreme_value.stan

/* 
Fit a GEV distribution in stan
Code is largely based on previous code by
Cameron Bracken: http://bechtel.colorado.edu/~bracken/tutorials/stan/
and Yenan Wu: http://discourse.mc-stan.org/t/troubles-in-rejecting-initial-value/1827
*/
functions{
  real gev_lpdf(vector y, real mu, real sigma, real xi) {
    vector[rows(y)] z;
    vector[rows(y)] lp;
    int N;
    N = rows(y);
    for(n in 1:N){
      z[n] = 1 + (y[n] - mu) * xi / sigma;
      lp[n] = (1 + (1 / xi)) * log(z[n]) + pow(z[n], -1/xi);
    }
    return -sum(lp) - N * log(sigma);
  }
}
data {
  int<lower=0> N;
  vector[N] y;
}

transformed data {
  real min_y;
  real max_y;
  real sd_y;
  min_y = min(y);
  max_y = max(y);
  sd_y = sd(y);
}
parameters {
  real<lower=-0.5,upper=0.5> xi;
  real<lower=0> sigma;
  // location has upper/lower bounds depending on the value of xi
  real<lower=if_else( xi > 0, min_y, negative_infinity()),
       upper=if_else( xi > 0, positive_infinity(), max_y )> mu;
}
model {
  // Priors -- these can be modified depending on your context
  real xi_plus_half;
  xi_plus_half = xi + 0.5;
  xi_plus_half ~ beta(1, 1);
  sigma ~ normal(sd_y, 10);
  mu ~ normal(0, 25);
  // Data Model
  y ~ gev(mu, sigma, xi);
}