In my Behavioral Economics class, one of the first things we learn about is risk preferences. Before teaching about this topic, I run a short, in-class experiment using the (free) Veconlab platform. This year, I got them to play several instances of the investment game, which is based on Gneezy and Potters (1997). Since a lot of my Behavioral Economics students are also my Econometrics students, I like to tie the two fields together as much as possible in this class. It is my research area, after all. So after yesterday’s class, I cleaned up the data and posted it with some questions as a problem set, and included in it for the Masters students the task of estimating a utility function from their own data. While I will focus on nonlinear least squares for their problem set (because I know they will be somewhat familiar with that), I couldn’t help busting out Stan and giving them the full Bayesian treatment as well. So what follows is a bit of a how-to guide for that.
The task
Students played ten rounds of the Gneezy and Potters (1997) task in Veconlab. In this task, a participant has to allocate an endowment \(m\) between a safe and a risky asset. The safe asset has a known, deterministic return, but the risky asset pays out a high amount with probability \(p\), and a low amount (less than the safe asset’s return) with probability \(1-p\). For yesterday’s experiment, I set the following parameters as constant:
- The endowment was set to \(m=0.25\).1
- The probability of the risky asset having the high return was \(p=0.5\).
- The low return of the safe asset was \(0.2\).
- The return of the safe asset was \(1\).
Then over the ten rounds I systematically varied the high return of the risky asset, call it \(\chi\in\{2.0, 2.4, 2.6, 2.8, \ldots , 4.0\}\). One of these rounds was randomly chosen for “payment” in extra credit points.
Some predictions
Here I will assume that my students each have constant relative risk-aversion (CRRA) expected utility preferences, which means that we can represent their expected utility from investing \(x\in[0,m]\) in the risky asset as:
\[ \begin{aligned} U_i(x)&=0.5 u_i(m-x+0.2x)+0.5 u_i(m-x+\chi x), \quad u_i(w)=\frac{w^{1-r_i}}{1-r_i} \end{aligned} \]
Although we won’t need it for estimation, it is informative to get the point prediction of this model by taking the first-order condition.
\[ \begin{aligned} 0&=-0.8\times0.5(m-x+0.2x)^{-r_i}+(\chi-1)0.5(m-x+\chi x)^{-r_i}\\ \frac{0.8}{\chi-1}&=\left(\frac{m+(\chi-1)x}{m-0.8x}\right)^{-r_i}\\ g_i(\chi)\doteq\left(\frac{\chi-1}{0.8}\right)^{\frac{1}{r_i}}&=\frac{m+(\chi-1)x}{m-0.8x}\\ mg_i(\chi)-0.8g_i(\chi)x&=m+(\chi-1)x\\ x&=m\frac{g_i(\chi)-1}{\chi-1+0.8g_i(\chi)} \end{aligned} \]
(provided that the fraction is in the unit interval, otherwise we have a corner solution).
So the predictions look like this:
library(tidyverse)
d<-expand_grid(
chi =c(2,3,4),
r = seq(-0.25,3, by=0.01)
) |>
mutate(
g = ((chi-1)/0.8)^(1/r),
x = 0.25*(g-1)/(chi-1+0.8*g)
) |>
mutate(
chi = paste("\u03c7 =",chi),
# where the FOC doesn't identify the maximum
x = ifelse(x<0,0.25,x),
x = ifelse(x>0.25, 0.25,x)
)
(
ggplot(d,aes(x=r,y=x,color=chi))
+geom_path(linewidth=1)
+theme_bw()
+xlab(expression(r[i]))
+scale_color_discrete(name="")
+ylab("optimal investment in risky asset")
)
So the task (with the parameters I chose) isn’t too good at identifying your \(r_i\) if you are risk-loving (\(r_i<0\)) or only a little bit risk-averse, but for a wide range of risk preferences in the risk-averse region the the model makes nice, interior solution predictions.
Or, with \(\chi\) on the horizontal axis, we could show the predictions for what people with different risk preferences would do:
d<-expand_grid(
chi = seq(2,4,by=0.01),
r = c(0, 0.4, 0.8, 1.2, 1.6, 2.0)
) |>
mutate(
g = ((chi-1)/0.8)^(1/r),
x = 0.25*(g-1)/(chi-1+0.8*g),
x = ifelse(x>=0.25, 0.25, x),
x = ifelse(r==0, 0.25, x),
r = paste("r =",r)
)
(
ggplot(d,aes(x=chi,y=x,color=as.factor(r)))
+geom_line(linewidth=1)
+ylab("optimal investment in risky asset")
+xlab(expression(chi~"(high return of risky asset)"))
+theme_bw()
+scale_color_discrete(name="")
+ylim(c(0,0.25))
)
A Stan model
Alongside the CRRA expected utility specification, I also assume a logistic choice rule. While students could in principle choose investments that weren’t increments of 0.01, fortunately they didn’t, so I assume that they noisily maximize over the discretized decision space \(x\in\{0.00, 0.01, 0.02, \ldots, 0.25\}\) with choice precision parameter \(\lambda_i\). Therefore I have two parameters to estimate for each participant: CRRA risk-aversion \(r_i\) and choice precision \(\lambda_i\).
Bayesian estimation requires priors for these parameters, which I set to:
\[ \begin{aligned} r_i&\sim N(0.27, 0.36^2)\\ \lambda_i&\sim \mathrm{LogNormal}(0,1) \end{aligned} \]
Here the prior for \(r_i\) I pulled directly from Bland (2025) (I’ve thought about this one a lot). This was calibrated to the distribution of choices in one of the treatments of Holt and Laury (2002). I have no idea about choice precision, so I just went with the prior for \(\lambda_i\) without much thought.2
Here is the Stan program I came up with to estimate my students’ risk preferences. Aside from the code to actually estimate the model, you can see that I have put a prediction in the generated quantities block. Here, we will compare (within-sample) the prediction made by the estimated model to the actual choices.
data {
int<lower=0> N; // number of observations
array[N] int<lower=0, upper=25> y; // integer amount invested
vector[N] Return; // high return of risky asset (chi in blog post)
vector[2] prior_r; // prior for r (normal, mean and sd)
vector[2] prior_lambda; // prior for lambda (lognormal, mean and sd)
}
transformed data {
real low_return = 0.2; // low return of risky asset
// the 26 numbers participants could choose between
vector[26] y_grid;
for (yy in 1:26) {
y_grid[yy] = yy-1;
}
}
parameters {
real r; // CRRA risk aversion
real<lower=0> lambda; // logit choice precision
}
model {
// prior
r ~ normal(prior_r[1],prior_r[2]);
lambda ~ lognormal(prior_lambda[1],prior_lambda[2]);
// likelihood
for (ii in 1:N) { // loop over each round of the task
// utility for each of the 26 possible choices
vector[26] U = 0.5*pow(25.0-y_grid+low_return*y_grid,1.0-r)/(1.0-r)+0.5*pow(25.0-y_grid+Return[ii]*y_grid,1.0-r)/(1.0-r);
target += categorical_logit_lpmf(y[ii]+1 | lambda*U);
}
}
generated quantities {
vector[N] prediction;
for (ii in 1:N) {
// utility for each of the 26 possible choices
vector[26] U = 0.5*pow(25.0-y_grid+low_return*y_grid,1.0-r)/(1.0-r)+0.5*pow(25.0-y_grid+Return[ii]*y_grid,1.0-r)/(1.0-r);
// choice probabilities for each of the 26 possible choices
vector[26] p = softmax(lambda*U);
// mean prediction.
prediction[ii] = p'*y_grid;
}
}And here is the R script I used to estimate the models. First I estimate a pooled model, then I estimate the model once for each student.
library(tidyverse)
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
GP<-"_posts/2026-01-21-estimating-my-students-risk-preferences/2026GPtask.csv" |>
read.csv() |>
mutate(
y = invested*100
)
model<-"_posts/2026-01-21-estimating-my-students-risk-preferences/GPindividual.stan" |>
stan_model()
D<-GP
d<-list(
N = dim(D)[1],
y = D$y |> as.integer(),
Return = D$Return,
prior_r = c(0.27,0.36),
prior_lambda = c(0,1)
)
Fit<-model |>
sampling(
data=d,
seed=42
)
ESTIMATES<-summary(Fit)$summary |>
data.frame() |>
rownames_to_column(var = "parameter") |>
mutate(
slice = "All data"
)
for (ii in unique(GP$ID)) {
print(ii)
D<-GP |> filter(ID==ii)
d<-list(
N = dim(D)[1],
y = D$y |> as.integer(),
Return = D$Return,
prior_r = c(log(0.5),1),
prior_lambda = c(0,1)
)
Fit<-model |>
sampling(
data=d,
seed=42,
# There were a few warnings, so I had to tune the sampler a bit
control = list(adapt_delta=0.99999),
iter = 5000,
chains = 8
)
ESTIMATES<-rbind(ESTIMATES,
summary(Fit)$summary |>
data.frame() |>
rownames_to_column(var = "parameter") |>
mutate(
slice = paste("ID =",ii)
)
)
}
ESTIMATES |>
write.csv("_posts/2026-01-21-estimating-my-students-risk-preferences/estimates.csv")
(
ggplot(ESTIMATES |> filter(parameter=="r"))
+theme_bw()
+geom_point(aes(x=mean,y=slice))
+geom_errorbar(aes(y=slice,xmin = X2.5., xmax = X97.5.))
+geom_vline(xintercept = 0,color="red",linetype="dashed")
+xlab("r")
+ylab("")
)Looking at the estimates
First, we can take a look at the posterior estimates of \(r_i\). Here I show the posterior means (dots) and 95% Bayesian credible regions (error bars, 2.5th-97.5th percentiles). We can see that there was a range of risk preferences, and unsurprisingly (because they only made ten decisions) there is a lot of posterior uncertainty in these estimates. The two credible regions that only span the risk-loving space represent two students who went all in for all ten decisions!
estimates<-"estimates.csv" |> read.csv()
data<-"2026GPtask.csv" |> read.csv()
(
ggplot(estimates |> filter(parameter=="r"))
+theme_bw()
+geom_point(aes(x=mean,y=slice))
+geom_errorbar(aes(y=slice,xmin = X2.5., xmax = X97.5.))
+geom_vline(xintercept = 0,color="red",linetype="dashed")
+xlab("r")
+ylab("")
)
Next, we can look at the pooled model’s predictions, and how they stack up against the average choices in the experiment:
chi<-c(2, 2.4, 2.6, 2.8, 3, 3.2, 3.4, 3.6, 3.8, 4)
data.summary<- data |>
group_by(Return) |>
summarize(
invested = mean(invested)
)
dplt<-estimates |>
filter(slice == "All data") |>
filter(grepl("prediction",parameter)) |>
mutate(
prediction.id = parse_number(parameter)
) |>
filter(
prediction.id<=10
) |>
mutate(
Return = chi[prediction.id]
) |>
full_join(
data.summary,
by = "Return"
) |>
mutate(
invested = invested*100
)
(
ggplot(dplt,aes(x=Return))
+geom_point(aes(y=invested))
+geom_ribbon(aes(ymin = X25., ymax = X75.),alpha=0.5, fill="red")
+theme_bw()
+ylim(c(0,25))
)
So the model is actually not too bad at organizing the aggregate data. But risk preferences are a model for the individual, not the aggregate! So we really should be looking at the individual choices and predictions from models estimated at the individual level.
r.estimates<-estimates |>
filter(parameter=="r" & slice !="All data") |>
mutate(
ID = parse_number(slice)
) |>
select(
ID, mean
) |>
rename(
r = mean
)
lambda.estimates<-estimates |>
filter(parameter=="lambda" & slice !="All data") |>
mutate(
ID = parse_number(slice)
) |>
select(
ID, mean
) |>
rename(
lambda = mean
)
dplt<-estimates |>
filter(slice != "All data") |>
filter(grepl("prediction",parameter)) |>
mutate(
ID = parse_number(slice),
prediction.id = parse_number(parameter),
Return = chi[prediction.id]
) |>
full_join(
data,
by = c("ID","Return")
) |>
mutate(
invested = invested*100
) |>
full_join(
r.estimates, by = "ID"
) |>
full_join(lambda.estimates, by="ID")|>
mutate(
facet.text = paste(ID, " r =",sprintf("%.2f",r), " \u03bb =",sprintf("%.2f",lambda))
)
(
ggplot(dplt,aes(x=Return))
+geom_point(aes(y=invested))
+facet_wrap(~facet.text)
+geom_ribbon(aes(ymin=X25., ymax = X75.), fill="red", alpha=0.5)
+geom_line(aes(y=mean),color="red")
+theme_bw()
+ylim(c(0,25))
)
Here we can see that the model fits almost perfectly for the students who went all in (IDs 4 and 5), and does fairly well for ID 6. For the others, things did not go so well. For me, this was somewhat disappointing for the predictions for IDs 1 and 3, as these students chose monotonically. It must be that it is hard for the CRRA model to match these predictions.
So can you do this with your class?
Absolutely! But if I were to do this again, I would either pick a different risk preference elicitation task, or have them make more than ten decisions. The estimates are quite noisy, and the model does not do a very good job at making predictions that look like what the monotonic participants are doing.
As for the feasibility of doing this, all of this ran on my laptop in about 10 minutes, and it all runs on free software (Veconlab, R, and Stan). I’ve given you all the code you need to do the estimation in this blog post.
If you have any questions about running this in your own clasroom, please just send me an email. Finally, if you’re interested in estimating structural models using Bayesian techniques, please check out my free online book!