-5

I'm trying to figure out why this doesn't work:

    data=read.csv("data_risk.csv")
    pa1 = c(data$pa1)
    pa2 = c(data$pa2)
    pb1 = c(data$pb1)
    pb2 = c(data$pb2)
    a1 = c(data$a1)
    a2= c(data$a2)
    b1 = c(data$b1)
    b2 = c(data$b2)
    yy=c(data$choice)

   crra=function(x,r){
    u=x^(1-r)/(1-r)
     return(u)
    }

   eua = c(pa1*crra(a1,r)+pa2*crra(a2,r))
   eub = c(pb1*crra(b1,r)+pb2*crra(b2,r))

   LL_all = c()

 R<-seq(0,1,0.01)
  for (r in R){
 eua = c(pa1*crra(a1,r)+pa2*crra(a2,r))
eub = c(pb1*crra(b1,r)+pb2*crra(b2,r))
 probA = eua/(eua+eub)
 total = ifelse(yy==1, probA, 1-probA)
  LL=log(prod(total))
 LL_all=c(LL_all,LL)
}

Right now every time I try and run it it says object r not found or error object R not found it works without the for loop just fine but when I add the for loop it all breaks down.

I'm trying to find the value of r that maximises someones utility given two choices. A decision maker chooses option A over B with probability a EUA/(EUA+EUB). In this example r is the risk aversion coefficient and x is the outcome of the lottery.

  1. pa1 = probability of event a1 happening
  2. pa2 = probability of event a2 happening
  3. pb1 = probability of event b1 happening
  4. pb2 = probability of event b2 happening
  5. a1,a2,b1,b2 = outcomes of events
  6. yy= indicator function that takes the value of 1 if lottery a is chosen and 0 otherwise

dataset:

 : task  pa1 a1  pa2 a2  pb1  b1  pb2 b2 choice
   1      0.34 24 0.66 59 0.42  47 0.58 64      0
   2      0.88 79 0.12 82 0.20  57 0.80 94      0
   3      0.74 62 0.26  0 0.44  23 0.56 31      1
   4      0.05 56 0.95 72 0.95  68 0.05 95      1
   5      0.25 84 0.75 43 0.43   7 0.57 97      0
   6      0.28  7 0.72 74 0.71  55 0.29 63      0
   7      0.09 56 0.91 19 0.76  13 0.24 90      0
   8      0.63 41 0.37 18 0.98  56 0.02  8      0
   9      0.88 72 0.12 29 0.39  67 0.61 63      1
   10    0.61 37 0.39 50 0.60   6 0.40 45      1
   11    0.08 54 0.92 31 0.15  44 0.85 29      1
   12    0.92 63 0.08  5 0.63  43 0.37 53      1
   13    0.78 32 0.22 99 0.32  39 0.68 56      0
   14    0.16 66 0.84 23 0.79  15 0.21 29      1
   15    0.12 52 0.88 73 0.98  92 0.02 19      0
   16    0.29 88 0.71 78 0.29  53 0.71 91      1
   17    0.31 39 0.69 51 0.84  16 0.16 91      1
   18    0.17 70 0.83 65 0.35 100 0.65 50      0
   19    0.91 80 0.09 19 0.64  37 0.36 65      1
   20    0.09 83 0.91 67 0.48  77 0.52  6      1
   21    0.44 14 0.56 72 0.21   9 0.79 31      1
   22    0.68 41 0.32 65 0.85 100 0.15  2      0
   23    0.38 40 0.62 55 0.14  26 0.86 96      0
   24    0.62  1 0.38 83 0.41  37 0.59 24      1
   25    0.49 15 0.51 50 0.94  64 0.06 14      0
   26    0.10 40 0.90 32 0.10  77 0.90  2      1
   27    0.20 40 0.80 32 0.20  77 0.80  2      1
   28    0.30 40 0.70 32 0.30  77 0.70  2      1
   29    0.40 40 0.60 32 0.40  77 0.60  2      1
   30    0.50 40 0.50 32 0.50  77 0.50  2      0
   31    0.60 40 0.40 32 0.60  77 0.40  2      0
   32    0.70 40 0.30 32 0.70  77 0.30  2      0
   33    0.80 40 0.20 32 0.80  77 0.20  2      0
   34    0.90 40 0.10 32 0.90  77 0.10  2      0
   35    1.00 40 0.00 32 1.00  77 0.00  2      0
Graham Thomson
  • 13
  • 2
  • 1
    [Provide minimal, reproducible, representative example(s) along with the desired end result. Use dput() for data and specify all non-base packages.](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) – Roman Dec 04 '18 at 23:01
  • 2
    So, screenshots are typically going to get a rather poor reaction from people on this site. They're really not very helpful in providing us the necessary information to help you. That said, I'm struggling to see any direct evidence in that screenshot that the object "r" is actually used successfully in a subsequent line. – joran Dec 04 '18 at 23:03

1 Answers1

0

The problem in the peace of code below after your definition of crra function:

eua = c(pa1*crra(a1,r)+pa2*crra(a2,r))
eub = c(pb1*crra(b1,r)+pb2*crra(b2,r))

Basically you are trying to use r variable before it's defined moreover it is a duplicate of the code inside the for-loop. If you comment out these two lines everything goes OK. Please see the code below:

data=read.table(text = " task  pa1 a1  pa2 a2  pb1  b1  pb2 b2 choice
   1      0.34 24 0.66 59 0.42  47 0.58 64      0
                2      0.88 79 0.12 82 0.20  57 0.80 94      0
                3      0.74 62 0.26  0 0.44  23 0.56 31      1
                4      0.05 56 0.95 72 0.95  68 0.05 95      1
                5      0.25 84 0.75 43 0.43   7 0.57 97      0
                6      0.28  7 0.72 74 0.71  55 0.29 63      0
                7      0.09 56 0.91 19 0.76  13 0.24 90      0
                8      0.63 41 0.37 18 0.98  56 0.02  8      0
                9      0.88 72 0.12 29 0.39  67 0.61 63      1
                10    0.61 37 0.39 50 0.60   6 0.40 45      1
                11    0.08 54 0.92 31 0.15  44 0.85 29      1
                12    0.92 63 0.08  5 0.63  43 0.37 53      1
                13    0.78 32 0.22 99 0.32  39 0.68 56      0
                14    0.16 66 0.84 23 0.79  15 0.21 29      1
                15    0.12 52 0.88 73 0.98  92 0.02 19      0
                16    0.29 88 0.71 78 0.29  53 0.71 91      1
                17    0.31 39 0.69 51 0.84  16 0.16 91      1
                18    0.17 70 0.83 65 0.35 100 0.65 50      0
                19    0.91 80 0.09 19 0.64  37 0.36 65      1
                20    0.09 83 0.91 67 0.48  77 0.52  6      1
                21    0.44 14 0.56 72 0.21   9 0.79 31      1
                22    0.68 41 0.32 65 0.85 100 0.15  2      0
                23    0.38 40 0.62 55 0.14  26 0.86 96      0
                24    0.62  1 0.38 83 0.41  37 0.59 24      1
                25    0.49 15 0.51 50 0.94  64 0.06 14      0
                26    0.10 40 0.90 32 0.10  77 0.90  2      1
                27    0.20 40 0.80 32 0.20  77 0.80  2      1
                28    0.30 40 0.70 32 0.30  77 0.70  2      1
                29    0.40 40 0.60 32 0.40  77 0.60  2      1
                30    0.50 40 0.50 32 0.50  77 0.50  2      0
                31    0.60 40 0.40 32 0.60  77 0.40  2      0
                32    0.70 40 0.30 32 0.70  77 0.30  2      0
                33    0.80 40 0.20 32 0.80  77 0.20  2      0
                34    0.90 40 0.10 32 0.90  77 0.10  2      0
                35    1.00 40 0.00 32 1.00  77 0.00  2      0", header = TRUE)

pa1 = c(data$pa1)
pa2 = c(data$pa2)
pb1 = c(data$pb1)
pb2 = c(data$pb2)
a1 = c(data$a1)
a2= c(data$a2)
b1 = c(data$b1)
b2 = c(data$b2)
yy=c(data$choice)

crra=function(x,r){
  u=x^(1-r)/(1-r)
  return(u)
}

# eua = c(pa1*crra(a1,r)+pa2*crra(a2,r))
# eub = c(pb1*crra(b1,r)+pb2*crra(b2,r))

LL_all = c()

R<-seq(0,1,0.01)
for (r in R){
  eua = c(pa1*crra(a1,r)+pa2*crra(a2,r))
  eub = c(pb1*crra(b1,r)+pb2*crra(b2,r))
  probA = eua/(eua+eub)
  total = ifelse(yy==1, probA, 1-probA)
  LL=log(prod(total))
  LL_all=c(LL_all,LL)
}
head(LL_all)

Output:

[1] -18.93759 -18.97863 -19.02000 -19.06170 -19.10374 -19.14611
Artem
  • 3,304
  • 3
  • 18
  • 41