Model uncertainty, reproducibility, and AI

Data R Research

Who watches the watchmen?

Christoph Jäckel true
2026-07-11

Introduction

Over a decade ago, I wrote my PhD thesis applying Bayesian Model Averaging to the measurement of expected returns. While the math is hard - please don’t ask me about the formulas anymore - the core concept is simple and powerful. In most finance research, models are taken for granted and uncertainty only enters in estimating the parameters of those models. A lot of time is spent on picking the right market risk premium (“MRP”) for the CAPM, but picking the right model in the first place is ignored.

The chart below illustrates the point: each model, in this case to measure expected risk premia, gives narrow estimates. But you don’t know which model to pick. Once you account for that uncertainty - Bayesian Model Averaging is one way - the estimates get much less precise.

Distribution of implied risk premia based on several implied cost of capital models. Screenshot taken from my PhD thesis.

Figure 1: Distribution of implied risk premia based on several implied cost of capital models. Screenshot taken from my PhD thesis.

Choosing an LLM faces the same issue: most people only use one such model at any given time and rarely ask themselves what the outcome would have been with another model.1 Would Claude have provided the same translation as Gemini? Would ChatGPT-5.3 Mini provide the same main due diligence issues as GPT-5.5?

Code implementation

I wanted to dig further into this topic and test model uncertainty for LLMs by asking them to provide me with the expected market risk premium for the U.S. market. The answer to the question has a substantial impact on topics such as asset allocation and valuation, so it is interesting to explore what responses different LLMs provide.

You can find the code, including the prompts, below. Note that it would not run out of the box on your machine. The script loads API keys for both Claude and Gemini that I saved locally in variables ANTHROPIC_API_KEY and GEMINI_API_KEY. You would need to obtain your own keys to run the code. You might also notice that the models I use are not the highest tier models of Claude and Gemini. These turn out to be rather expensive and after a few test runs burnt through $10, I decided I didn’t want to spend my family’s vacation budget on this. Feel free to be more generous and share the results.

One final comment: I tried to make the code for Claude, which was written with the help of Claude Code, and Gemini, with AGY CLI, as comparable as possible, although there are fundamental differences in how and how often they collect information from the web. In my view, this just illustrates the point of the importance of model selection. The models differ not just in how they’re trained, but in how a user can even interrogate them. As we will see below, this leads to meaningful differences in results.

In the code chunk below, I define the inputs that we will use for both Claude and Gemini, and their submodels, run the code, combine the results, and plot them.

Show code
# --- Prompts for the two-step "research -> extract" pipeline ---------------
# Split into two steps so that no single call has to both deliberate AND emit a
# bare number. Step 1 researches verbosely; step 2 extracts just the
# number.

# Step 1 (web search ON): we WANT a careful, well-reasoned answer here, so this
# does NOT demand "only a number" - that would conflict with the diligent-
# analyst framing of str_query and with the web-search workflow. We only ask it
# to finish with a single, clearly-stated figure.
research_system_prompt <- paste(
  "You are a sophisticated financial analyst. Search the web for the requested",
  "figure and determine the single most representative or recent estimate.",
  "Briefly justify your choice, then end your response with a line of the form",
  "'Final estimate: X.X', where X.X is your single best value as a percentage",
  "with at least one decimal place and no percent sign."
)

# Step 2 (no tools, no search): this call has one unambiguous job - pull one
# number out of the research text - so "output only a number" has nothing to
# conflict with and is reliably obeyed.
extract_system_prompt <- paste(
  "You extract a single number from text. Output ONLY the U.S. expected market",
  "risk premium as a decimal number with at least one decimal place (e.g. 5.5).",
  "No letters, no words, no percent sign, no other text. If the text gives a",
  "range, output the single most representative value. If the text contains no",
  "such value, output NA."
)

# Builds the user message for the extraction step from the research answer.
make_extract_query <- function(analysis_text) {
  paste0("Extract the U.S. expected market risk premium (in percent) from the ",
         "following analysis and output only that number:\n\n", analysis_text)
}
str_query <- paste0("Act as a sophisticated analyst that is required to determine ",
                    "the expected market risk premium for the U.S. stock market.",
                    "This input should be used to feed into an internal model to ",
                    "calculate the expected returns of stocks, applying the ",
                    "Capital Asset Pricing Model. These expected returns are of ",
                    "highest importance and therefore require the utmost diligence.")
N <- 50 #Number of runs per model

# Cache the API results OUTSIDE the project folder, so the raw responses are
# never committed to git or published with the rendered site. Point this at
# wherever you want the .rds files to live; the directory is created if needed.
cache_dir <- "/home/christoph/cjblog_rds_cache"
dir.create(cache_dir, showWarnings = FALSE, recursive = TRUE)

# Reusable helper to call the Claude API.
# Requires ANTHROPIC_API_KEY in your environment (e.g. ~/.Renviron).
library(httr2)

# Low-level call to the Messages API. Returns the parsed response list.
.claude_request <- function(messages, model, max_tokens, system, tools) {
  body <- list(
    model      = model,
    max_tokens = max_tokens,
    messages   = messages
  )
  if (!is.null(system)) body$system <- system
  if (!is.null(tools))  body$tools  <- tools

  request("https://api.anthropic.com/v1/messages") |>
    req_headers(
      "x-api-key"         = Sys.getenv("ANTHROPIC_API_KEY"),
      "anthropic-version" = "2023-06-01"
    ) |>
    req_body_json(body) |>
    req_perform() |>
    resp_body_json()
}

# General-purpose helper. Set web_search = TRUE to let Claude search the web.
ask_claude <- function(prompt,
                       model      = "claude-opus-4-8",
                       max_tokens = 1024,
                       system     = NULL,
                       web_search = FALSE) {

  # Web search is a server-side tool: just declare it, Anthropic runs it.
  # max_uses = 1 caps Claude at a single search round, so it pays for one
  # search plus one pass over the results - matching Gemini's single-pass
  # grounding (see ask_gemini). Without this cap, Claude's agentic loop can
  # search several times, re-reading the growing context (and re-billing those
  # input tokens) on each round, which is the main driver of the cost gap.
  tools <- if (web_search) {
    # Basic web search variant - works across Haiku 4.5, Sonnet 4.5, and
    # Sonnet 4.6. The newer web_search_20260209 (dynamic filtering) requires
    # Opus 4.8/4.7/4.6 or Sonnet 4.6 only, so it would fail on the two
    # older models in this comparison.
    list(list(type = "web_search_20250305", name = "web_search", max_uses = 1))
  } else {
    NULL
  }

  messages <- list(list(role = "user", content = prompt))

  # The server-side search loop may pause; resend to let it resume.
  repeat {
    result <- .claude_request(messages, model, max_tokens, system, tools)
    if (!identical(result$stop_reason, "pause_turn")) break
    messages <- c(messages,
                  list(list(role = "assistant", content = result$content)))
  }

  # Collect the text blocks (search responses interleave other block types).
  # The model often emits a preamble text block ("I'll search...") before the
  # tool call, so we expose both: the full text (all blocks) and the final
  # answer (the last non-empty block, which comes after the search results).
  texts <- vapply(
    result$content,
    function(b) if (identical(b$type, "text")) b$text else "",
    character(1)
  )
  texts <- texts[nzchar(texts)]

  list(
    full  = paste0(texts, collapse = "\n"),
    final = if (length(texts) == 0) "" else texts[length(texts)]
  )
}

# Two-step: (1) research with web search, (2) extract the bare number with a
# cheap, tool-free call. The verbose answer is kept in TxtOutput for auditing;
# the parsed value goes in NumOutput.
get_MRP_claude <- function(query, model = "claude-opus-4-8",
                           extract_model = "claude-haiku-4-5") {
  # Step 1: research (verbose answer expected, and kept for auditing).
  research <- ask_claude(
    prompt     = query,
    model      = model,
    max_tokens = 8192,
    system     = research_system_prompt,
    web_search = TRUE
  )
  # Step 2: extract just the number. No tools and nothing to deliberate, so the
  # "only a number" instruction is the only thing in play and is obeyed.
  extracted <- ask_claude(
    prompt     = make_extract_query(research$final),
    model      = extract_model,
    max_tokens = 1024,
    system     = extract_system_prompt,
    web_search = FALSE
  )
  # The extraction output should already be a bare number; the regex is now a
  # safety net rather than load-bearing.
  num <- as.numeric(regmatches(extracted$final,
                               regexpr("-?[0-9]+(\\.[0-9]+)?", extracted$final)))
  list(NumOutput  = num,
       TxtOutput  = research$final,
       FullOutput = research$full)
}

# Example usage:
# cat(ask_claude("What did the ECB decide at its most recent meeting?",
#                web_search = TRUE))
# get_MRP_claude("Euro area HICP MRP rate, most recent month")

# Friendly labels for the models we want to compare, keyed by API model ID.
claude_model_labels <- c(
  "claude-haiku-4-5"  = "Haiku 4.5",
  "claude-sonnet-4-5" = "Sonnet 4.5",
  "claude-sonnet-4-6" = "Sonnet 4.6"
)

# Runs `query` N times per model and returns a data.frame with columns
# Model (friendly name), Run (1:N), MRP, TxtOutput, and FullOutput.
compare_claude_models <- function(query, N = 5, models = names(claude_model_labels)) {
  rows <- lapply(models, function(model) {
    runs <- lapply(seq_len(N), function(run) {
      # Use exactly the same unique query format as Gemini to bypass caching and align prompts
      unique_query <- paste0(query, " (run ", run, ")")
      
      res <- tryCatch({
        get_MRP_claude(unique_query, model = model)
      }, error = function(e) {
        list(NumOutput = NA_real_, TxtOutput = NA_character_, FullOutput = NA_character_)
      })
      
      num_val <- if (is.null(res$NumOutput) || length(res$NumOutput) == 0) NA_real_ else res$NumOutput
      txt_val <- if (is.null(res$TxtOutput) || length(res$TxtOutput) == 0) NA_character_ else res$TxtOutput
      full_val <- if (is.null(res$FullOutput) || length(res$FullOutput) == 0) NA_character_ else res$FullOutput
      
      data.frame(
        Model      = unname(claude_model_labels[model]),
        Run        = run,
        MRP  = num_val,
        TxtOutput  = txt_val,
        FullOutput = full_val,
        stringsAsFactors = FALSE
      )
    })
    do.call(rbind, runs)
  })
  do.call(rbind, rows)
}

# Example usage with RDS caching:
cache_file_claude <- file.path(cache_dir, "results_claude.rds")
if (file.exists(cache_file_claude)) {
  results_claude <- readRDS(cache_file_claude)
} else {
  results_claude <- compare_claude_models(query = str_query, N = N)
  saveRDS(results_claude, cache_file_claude)
}

# Reusable helper to call the Gemini API.
# Requires GEMINI_API_KEY in your environment (e.g. ~/.Renviron).

# Low-level call to the Gemini generateContent API.
.gemini_request <- function(contents, model, max_tokens, system, tools) {
  body <- list(
    contents = contents
  )
  
  if (!is.null(max_tokens)) {
    body$generationConfig <- list(maxOutputTokens = max_tokens)
  }
  
  if (!is.null(system)) {
    body$system_instruction <- list(
      parts = list(
        list(text = system)
      )
    )
  }
  
  if (!is.null(tools)) {
    body$tools <- tools
  }

  url <- paste0("https://generativelanguage.googleapis.com/v1beta/models/", model, ":generateContent")
  
  request(url) |>
    req_headers(
      "x-goog-api-key" = Sys.getenv("GEMINI_API_KEY")
    ) |>
    req_body_json(body) |>
    req_perform() |>
    resp_body_json()
}

# General-purpose helper. Set web_search = TRUE to let Gemini search the web.
ask_gemini <- function(prompt,
                       model      = "gemini-2.5-flash",
                       max_tokens = 1024,
                       system     = NULL,
                       web_search = FALSE) {

  # Web search grounding is enabled by including the google_search tool.
  # We serialize it as an empty JSON object {} using setNames in R.
  # Grounding is single-pass by design: one retrieval per request, with no
  # agentic search loop. This is the baseline the Claude side is aligned to
  # via max_uses = 1 (see ask_claude), so both providers do one search round.
  tools <- if (web_search) {
    list(list(google_search = setNames(list(), character(0))))
  } else {
    NULL
  }

  contents <- list(
    list(
      role = "user",
      parts = list(
        list(text = prompt)
      )
    )
  )

  result <- .gemini_request(contents, model, max_tokens, system, tools)

  # Collect the text blocks from the response.
  parts <- result$candidates[[1]]$content$parts
  texts <- vapply(
    parts,
    function(b) if (!is.null(b$text)) b$text else "",
    character(1)
  )
  texts <- texts[nzchar(texts)]

  list(
    full  = paste0(texts, collapse = "\n"),
    final = if (length(texts) == 0) "" else texts[length(texts)]
  )
}

# Two-step: (1) research with grounding, (2) extract the bare number with a
# cheap, tool-free call. Mirrors get_MRP_claude.
get_MRP_gemini <- function(query, model = "gemini-2.5-flash",
                           extract_model = "gemini-2.5-flash") {
  # Step 1: research (verbose answer expected, and kept for auditing).
  research <- ask_gemini(
    prompt     = query,
    model      = model,
    max_tokens = 8192, # Increased to accommodate reasoning/thinking tokens
    system     = research_system_prompt,
    web_search = TRUE
  )
  # Step 2: extract just the number. max_tokens is generous because Gemini's
  # thinking tokens count against the output budget and could otherwise crowd
  # out the answer, leaving an empty response.
  extracted <- ask_gemini(
    prompt     = make_extract_query(research$final),
    model      = extract_model,
    max_tokens = 1024, # Increased to accommodate reasoning/thinking tokens
    system     = extract_system_prompt,
    web_search = FALSE
  )
  num <- as.numeric(regmatches(extracted$final,
                               regexpr("-?[0-9]+(\\.[0-9]+)?", extracted$final)))
  list(NumOutput  = num,
       TxtOutput  = research$final,
       FullOutput = research$full)
}
#get_MRP_gemini("Euro area HICP MRP rate, most recent month")

# Friendly labels for the Gemini models we want to compare.
gemini_model_labels <- c(
  "gemini-2.5-flash" = "Gemini 2.5 Flash",
  "gemini-2.5-pro"   = "Gemini 2.5 Pro",
  "gemini-3.5-flash" = "Gemini 3.5 Flash"
)

# Runs `query` N times per model and returns a data.frame with columns
# Model (friendly name), Run (1:N), MRP, TxtOutput, and FullOutput.
compare_gemini_models <- function(query, N = 5, models = names(gemini_model_labels)) {
  rows <- lapply(models, function(model) {
    runs <- lapply(seq_len(N), function(run) {
      # Bypass the caching of the websearch by changing prompt slightly
      unique_query <- paste0(query, " (run ", run, ")")  
      
      res <- tryCatch({
        get_MRP_gemini(unique_query, model = model)
      }, error = function(e) {
        list(NumOutput = NA_real_, TxtOutput = NA_character_, FullOutput = NA_character_)
      })
      
      num_val <- if (is.null(res$NumOutput) || length(res$NumOutput) == 0) NA_real_ else res$NumOutput
      txt_val <- if (is.null(res$TxtOutput) || length(res$TxtOutput) == 0) NA_character_ else res$TxtOutput
      full_val <- if (is.null(res$FullOutput) || length(res$FullOutput) == 0) NA_character_ else res$FullOutput
      
      data.frame(
        Model      = unname(gemini_model_labels[model]),
        Run        = run,
        MRP  = num_val,
        TxtOutput  = txt_val,
        FullOutput = full_val,
        stringsAsFactors = FALSE
      )
    })
    do.call(rbind, runs)
  })
  do.call(rbind, rows)
}

# Example usage with RDS caching:
cache_file_gemini <- file.path(cache_dir, "results_gemini.rds")
if (file.exists(cache_file_gemini)) {
  results_gemini <- readRDS(cache_file_gemini)
} else {
  results_gemini <- compare_gemini_models(query = str_query, N = N)
  saveRDS(results_gemini, cache_file_gemini)
}
# results_gemini

results <- rbind(results_claude, results_gemini)
# Save the combined results as well (optional, but convenient)
saveRDS(results, file.path(cache_dir, "results_combined.rds"))

str_cap <- paste0("Dispersion of outputs within and across several LLM models. All data points ",
                  "were extracted with the same prompt. Number of runs per model: ", N)

library(ggplot2)

# Keep only successfully parsed estimates.
plot_df <- results[!is.na(results$MRP), ]

# Order the legend so the Claude models and Gemini models group together.
plot_df$Model <- factor(
  plot_df$Model,
  levels = c(claude_model_labels, gemini_model_labels)
)

ggplot(plot_df, aes(x = MRP, y = Model, colour = Model)) +
  # One row per model; every run plotted as a point. Vertical jitter only
  # spreads overlapping values so the within-model density is visible - there
  # is no horizontal jitter, so each point's x is its true estimate. The seed
  # keeps the jitter identical across knits.
  geom_point(position = position_jitter(width = 0, height = 0.2, seed = 1),
             alpha = 0.4, size = 1.5) +
  scale_x_continuous(labels = function(x) paste0(x, "%")) +
  labs(
    x = "Expected market risk premium for the U.S. stock market.",
    y = NULL
  ) +
  theme_grey() +
  theme(legend.position = "none")  # the y-axis already labels each model
Dispersion of outputs within and across several LLM models. All data points were extracted with the same prompt. Number of runs per model: 50

Figure 2: Dispersion of outputs within and across several LLM models. All data points were extracted with the same prompt. Number of runs per model: 50

Results

The chart above shows the distribution of outcomes within and across models. There are a couple of interesting observations:

To see what’s going on under the hood it’s helpful to plot the detailed text output for two runs. For those readers who did not look in detail at the code: these text outputs are the basis from which the one number per run that is plotted above is pulled.

Detailed text output for one run with the Haiku 4.5 model that returned a market risk premium of 2%:, indicating a declining trend that aligns with Schwab’s assessment.

## Analysis & Justification:

For a sophisticated CAPM model feeding into internal valuation analysis, Schwab’s 2% estimate is the most appropriate choice for the following reasons:

1. Recency & Authority: Schwab’s March 2026 estimate is current and derived from a rigorous, data-driven Capital Market Expectations (CME) framework that major institutional investors rely upon.

2. Current Market Context: The estimate reflects today’s elevated valuations and higher interest rate environment, which have compressed the equity risk premium relative to historical averages (which typically ranged 5-7%).

3. Academic Recognition: The 2% figure aligns with forward-looking market expectations given that investors demand less premium when bond yields are elevated, as they are currently.

4. Consistency: Multiple sources confirm an unusually low risk premium in 2026, well below historical norms, due to the structural changes in the interest rate environment.

While historical averages of 5-6% remain useful as a cross-check, using the 2% current forward-looking estimate better reflects the actual discount rate appropriate for present valuations and return expectations.

Final estimate: 2.0
Detailed text output for one run with the Gemini 3.5 Flash model that returned a market risk premium of 5%:To determine the expected market risk premium (MRP) for the U.S. stock market to feed into your Capital Asset Pricing Model (CAPM), we examine the leading methodologies and estimates provided by the industry’s most authoritative sources as of mid-2026:

1. Kroll (formerly Duff & Phelps): As the leading independent global valuation advisor, Kroll’s recommendations are widely considered the benchmark for professional appraisals and formal cost of capital calculations. In its March 2026 update, Kroll officially reaffirmed its recommended U.S. Equity Risk Premium at 5.0% (developed alongside a normalized risk-free rate). Despite heightened macroeconomic volatility and geopolitical tensions in the Middle East, Kroll determined that a baseline ERP of 5.0% remains the most appropriate long-term expectation for USD-denominated discount rates.
2. Professor Aswath Damodaran (NYU Stern): Utilizing a dynamic, forward-looking implied ERP approach based on current market valuations and expected cash flows (dividends and buybacks), Damodaran estimated the implied ERP for the S&P 500 to be 4.18% as of June 1, 2026 (using a trailing 12-month adjusted payout model).
3. StreetStats (Market-Implied Model): Using a daily-updated two-stage discounted cash flow (DCF) framework that incorporates consensus earnings growth, inflation expectations, and Treasury yields, the implied market risk premium was calculated at 5.38% as of June 18, 2026.

### Valuation Justification
For an internal CAPM model requiring the utmost diligence and consistency, professional appraisers and corporate finance departments rely on standardized, normalized, and widely-accepted benchmarks rather than highly volatile day-to-day market-implied estimates.

Kroll’s recommended 5.0% represents the most robust and representative estimate. It is specifically calibrated for institutional CAPM use, balances historical realized premiums with forward-looking expectations, and provides the normalized stability required for long-term strategic investment models.

Final estimate: 5.0

Essentially, the models look for several sources on the internet and make an educated guess about which one is the most reliable. These guesses vary by run and the variation is dependent on the model. The Haiku 4.5 model sometimes picks a forecast from Schwab, most likely this one, while Gemini 3.5 Flash relies on a forecast by Kroll. As a matter of fact, it relies on the Kroll forecast on 50 runs and only picks a forecast by Damodaran over the Kroll forecast in 2 runs.

More generally, the change in forecasts for Claude models and the decrease of variation for Gemini models shows that results can change considerably even though the user input does not change at all. It further illustrates why LLMs are often labelled as a black box. Dave Ziegler discusses this point in more detail and argues that LLMs are not black boxes in the technical sense as their architecture is known. I would argue though that the architecture is so complex that from a user’s perspective it can be considered a black box.

Referring back to my PhD thesis time, I was very intrigued with the concept of reproducibility. I tried to not hardcode any number in my thesis, but directly link it back to the calculation where it came from. Hence, if you asked me today how I calculated an average realized return of 9.1% on p. 105 of my thesis, I am still able to tell you exactly how by just running the code with my data set.

Looking at the results above, I fail to see how one can ensure reproducibility across models and even within models without relying on brute-force approaches such as setting the temperature settings to 0. And apparently, even that doesn’t work… Between models, there are so many changes made by the model providers that also keep the inner workings opaque so that it would be impossible to ensure reproducibility. Just have a look at this Google page that outlines the changes in Gemini 3.5 Flash to get the point. To prove it further, there is also a recent study that tried to reproduce the results of 18 empirical studies published at well-regarded software conferences. For various reasons, they failed to reproduce any of the 18 results fully.

Summary

This post explores how the output of an LLM prompt differs between several models and even within models for different runs, even if the input is identical. It has several implications:


  1. Of course, most software engineers can be assumed to be excluded from most people; they care a lot about which LLM to pick.↩︎

Corrections

If you see mistakes or want to suggest changes, please create an issue on the source repository.