3

I'm trying to use R to estimate E[u(X)] where u is a utility function and X is a random variable. More specifically, I want to be able to rank E[u(X)] and E[u(Y)] for two random variables X and Y -- only the ranking matters.

My problem is that u(x) = -exp(-sigma * x) for some sigma > 0, and this converges very rapidly to zero. So I have many cases where I expect, say, E[u(X)] > E[u(Y)], but because they are so close to zero, my simulation cannot distinguish them.

Does anyone have any advice for me?

I am only interested in ranking the two expected utilities, so u(x) can be replaced by any u.tilde(x) = a * u(x) + b, where a > 0 and b can be any number.

Below is an example where X and Y are both normal (in which case I think there is a closed form solution, but pretend X and Y have complicated distributions that I can only simulate from).

get.u <- function(sigma=1) {
    stopifnot(sigma > 0)
    utility <- function(x) {
        return(-exp(-sigma * x))
    }
    return(utility)
}

u <- get.u(sigma=1)
curve(u, from=0, to=10)  # Converges very rapidly to zero

n <- 10^4
x <- rnorm(n, 10^4, sd=10)
y <- rnorm(n, 10^4, sd=10^3)
mean(u(x)) == mean(u(y))  # Returns True (they're both 0), but I expect E[u(x)] > E[u(y)]

## An example of replacing u with a*u + b
get.scaled.u <- function(sigma=1) {
    stopifnot(sigma > 0)  # Risk averse
    utility <- function(x) {
        return(-exp(-sigma * x + sigma * 10^4))
    }
    return(utility)
}

u <- get.scaled.u(sigma=1)
mean(u(x)) > mean(u(y))  # True as desired

x <- rnorm(n, 10^4, sd=10^3)
y <- rnorm(n, 10^4, sd=2*10^3)
mean(u(x)) > mean(u(y))  # False again -- they're both -Inf

Is finding a clever way to scale u the correct way to deal with this problem? For example, suppose X and Y both have bounded support -- if I know the bounds, how can I scale u to guarantee that a*u + b will be neither too close to -Inf, nor too close to zero?

Edit: I didn't know about multiple precision packages. Rmpfr is helpful:

library(Rmpfr)
x.precise <- mpfr(x, 100)
y.precise <- mpfr(y, 100)
mean(u(x.precise)) > mean(u(y.precise))  # True
Adrian
  • 3,138
  • 2
  • 28
  • 39
  • Some tips: first version: `u(-710) == -Inf` `u(746) == 0`. Scaled version: `u(10^4-710) == -Inf` `u(10^4+746) == 0`. – Ferdinand.kraft Dec 01 '13 at 23:16
  • Thank you @Ferdinand.kraft, that makes sense. So, in the case sigma=1 ie u(x) = -exp(-x), it seems that I can never deal with an X whose support has length >= 710+745 = 1455. Is that correct? – Adrian Dec 02 '13 at 21:05
  • And if so, should I look for something along the lines of http://stackoverflow.com/questions/9559346/deal-with-overflow-in-exp-using-numpy? – Adrian Dec 02 '13 at 21:05

0 Answers0