The answer should be 117.4, I'm getting 9982.3... Not sure what the problem is but here's my code:
def util(c,p,alpha):
mu = 0
for i in range(0,len(c)):
m = p[i]*(c[i]**(1-alpha))
mu += m
return mu**(1/(1-alpha))
omega = 0.495
c = np.array([100,200,1000])
p = np.array([omega, omega, 1-2*omega])
alpha = 5
EDIT: I'm not sure if there is an error with my math or with the function I wrote, I am asking if my math fits the code that I wrote.
I'm solving this equation for mu: U(mu) = E[U(c)] with payoffs c, and probability distribution p, as above. U(c) has the form c^(1-alpha)/(1-alpha).
U(mu) = mu^(1-alpha)/(1-alpha) = E[U(c)] = (omega*c1^(1-alpha)+omega*c2^(1-alpha)+(1-2*omega)*c3^(1-alpha))/(1-alpha)
=> mu = (omega*c1^(1-alpha)+omega*c2^(1-alpha)+(1-2*omega)*c3^(1-alpha))^(1/(1-alpha))
Your main problem is that Python is doing integer division. Python 2 does integer division with /, unless you do from __future__ import division (see PEP 238). So you need to change at least one of your operands to a floating-point value. You can do this by setting alpha = 5.0 instead of alpha = 5. Or you can write 1.0 - alpha instead of 1 - alpha.
Also, you can make your code a bit more compact by using numpy's vectorized operations. Your util function can be shortened to:
def util(c, p, alpha):
mu = np.sum(p * (c ** (1.0 - alpha)))
return mu ** (1 / (1.0 - alpha))