Case Study: Stochastic Production and Cost Frontiers

Modeling
1

I’m using Stan for some different projects and decided that it might be useful to write out the basic stochastic frontier models as a case study. These are “composed error models” and used for efficiency analysis in econometrics. The models are not that complicated but I would appreciate any feedback you can provide. In particular, is there a more elegant way of combining a user-defined function and using the loo package? I’ve found Stan to be a useful piece of software for my work, but I’m never quite sure that I’m using best practices.

Check the repository for the case study at
https://bitbucket.org/aframsey/stancase_frontier/src/master/
or view a PDF of the final output. The PDF is a little wonky.
spfstan.pdf (144.4 KB)

If there’s interest in this case study, I may update with a few of the more advanced models for panel data. Or hierarchical models. But that might generate a whole new set of questions for the Stan forums!

2 Likes
2

Just for concreteness, I’m pasting in one of the models from the linked repo:

functions {
   
   real hnormal_lpdf(vector y, real sigma) {

      int N = rows(y);

      vector[N] prob;
      real lprob;
      for (i in 1:N) {
         prob[i] = log(2*sigma)-log(pi())-(y[i]^2*sigma^2)/pi();
      }
      lprob = sum(prob);
      return lprob;
   }

   real nhnormal_pdf(real y, row_vector x, real alpha, vector beta, real sigma, real lambda) {
      real prob;
      real sigs = sqrt(sigma^2 + lambda);
      real delt = sqrt(lambda)/sigma;
      real epsi = y - alpha - x * beta;
      prob = (0.5) * log(2 / pi()) - log(sigs) + normal_lcdf(epsi*delt/sigs|0,1) - (epsi^2)/(2*sigs^2);
      return prob;
   }
}

data {
   int<lower=0> N; // Number of observations
   int<lower=0> K; // Number of factors of cost function
   matrix[N, K] x; // Matrix of factors of cost function
   vector[N] y;    // Vector of costs
}
parameters {
   real alpha;            // Intercept
   vector[K] beta;        // Coefficients of the cost function
   real<lower=0> sigma;   // Error standard deviation
   real<lower=0> lambda;  // Parameter of half-normal inefficiency distribution
   vector<lower=0>[N] u;  // Inefficiency terms
}
model {
   alpha ~ normal(0, 100);
   beta ~ normal(0, 100);
   sigma ~ inv_gamma(0.001, 0.001);
   lambda ~ gamma(1, 1/37.5);
   u ~ hnormal(lambda);
   y ~ normal(alpha + x * beta + u, sigma);
}
generated quantities {
   vector[N] eff;
   vector[N] log_lik;
   eff = exp(-u);
   for (i in 1:N) {
      log_lik[i] = nhnormal_pdf(y[i], x[i], alpha, beta, sigma, lambda);
   }
}

First, we don’t recommend those diffuse gamma priors—there’s a discussion in the Stan manual in the chapter on regression. We generally recommend at least weakly informative priors that indicate the scale of the parameters.

You can vectorize for efficiency and you don’t need all those intermediates. This,

   real hnormal_lpdf(vector y, real sigma) {

      int N = rows(y);

      vector[N] prob;
      real lprob;
      for (i in 1:N) {
         prob[i] = log(2*sigma)-log(pi())-(y[i]^2*sigma^2)/pi();
      }
      lprob = sum(prob);
      return lprob;
   }

can be reduced to a much more efficient one liner

real hnormal_lpdf(vector y, real sigma) {
  return rows(y) * (log(2 * sigma) - log(pi()))
     - dot_self(y) / pi() * square(sigma);
}

with only N + 3 multiplications, one division, and two subtractions. I took the additional measure of guessing that y might be data, so that the dot-self and division by pi() are together so as to minimize operations involving parameters (which is where the real cost is because of the need to store the expressions and compute derivatives).

2 Likes
3

Hi

I know that your post about Stochastic Productions frontiers is too old, but I’m interested in this kind of models. Are your project alive?