1

I try to estimate the relative rate of risk aversion gmma of a CRRA function using R. The likelihood maximisation either does not converge or produces an unrealistic estimate. In order to make the problem I simulated choice data. I suspect a problem in the way I generate the data but I cannot spot it.

rm(list = ls())

library(haven)
library(tidyverse)

#enter utility function
utility <- function(x, gmma){
  if(gmma==1){
    return(log(x))
  }else{
    return((x^(1-gmma)-1)/(1-gmma))
  }
}


#function to calculate negative log likelihood
nll<- function(choice, p, gmma, ce, loutc){
  #read me
  #choice is an nxk matrix where n is the number of observations and k is the number of decisions
  #p is a vector length number of options in each lottery
  #ce is the certainty eqivalent. to keep things simple
  #gmma is the relative rate of risk aversion which is estimated
  #utility is the utility function that has to be specified before
  #loutc is a kxm matrix where k is the number of decisions and m is the number of options. 
  
  
  #calculate expected utility of lottery
  for (i in c(1:ncol(loutc))) {
    if(i==1){
      EU_L<-p[i]*utility(loutc[,i], gmma = gmma)
    } 
    else{
      EU_L<-EU_L+p[i]*utility(loutc[,i], gmma = gmma)
    }
  }
  
  #calculate expected utility of certain choice
  EU_ce <-utility(ce, gmma = gmma)
  
  #get fechner index
  euDiff<-EU_L-EU_ce
  
  #evaluate likelihood using observed choices
  eEuDiff<-matrix(NA,nrow=nrow(choice), ncol=ncol(choice))
  #do the inner product
  for (i in c(1:ncol(choice))){
    eEuDiff[,i]<-as.vector(as.numeric((choice[,i])))*euDiff[i]
  }
  #calculate matrix of densities
  ll_mat<- apply(eEuDiff, c(1,2), dnorm, mean = 0, sd = 1, log = TRUE)
  #calculate log likelihood
  ell_vec=apply(ll_mat, 1, sum)
  
  
  print("nll")
  print(nll)
  
  nll<- -sum(ell_vec, na.rm=TRUE)
  
  return(nll)
}

#simulate choice
choice_funct<-function(n, loutc, ce, p, gmma){
  #calculate expected utility of lottery
  for (i in c(1:ncol(loutc))) {
    if(i==1){
      EU_L<-p[i]*utility(loutc[,i], gmma = gmma)
    } 
    else{
      EU_L<-EU_L+p[i]*utility(loutc[,i], gmma = gmma)
    }
  }
  #Utility of certain option
  EU_ce <-utility(ce, gmma = gmma)
  
  #get fechner index
  
  euDiff<-EU_L-EU_ce
  #create som noise
  sd_ep<-abs(.5*sd(euDiff))
  epsilon<-cbind(rep(1,n),
                 matrix(rnorm(n=(n*length(euDiff)), mean=0, sd=sd_ep), ncol = length(euDiff))
  )
  euDiffMat<-rbind(euDiff,diag(rep(1, length(euDiff))))
  #add noise to utility difference
  euDiffRand<-epsilon %*% euDiffMat
  print(euDiffRand)
  #create choice result matrix
  chsimul<-ifelse(euDiffRand >0, 1, -1)
  return(chsimul)
}


p<-c(.5,.5)
loutc<-cbind(rep(1000,5), c(1500,1800,2000,2300,3000))
ce=1600
chs<-choice_funct(n=30,loutc = loutc, ce=ce, p=p, gmma=.5)
chs




optim(par = .5,fn= nll, choice=chs, p=p, ce=ce, loutc = loutc, method = "Brent", lower = -3, upper = 5)

When I run this optimization the result is

<bytecode: 0x0000025510793540>
$par
[1] 3.521294

$value
[1] 137.8408

$counts
function gradient 
      NA       NA 

$convergence
[1] 0

$message
NULL

Wich is ages away from the value of .5 for gamma that I used when generating the data. There must be a problem in the lower part of the nll function. For clarification I repeat the part of code where I expec the error to be.

  #evaluate likelihood using observed choices
  eEuDiff<-matrix(NA,nrow=nrow(choice), ncol=ncol(choice))
  #do the inner product
  for (i in c(1:ncol(choice))){
    eEuDiff[,i]<-as.vector(as.numeric((choice[,i])))*euDiff[i]
  }
  #calculate matrix of densities
  ll_mat<- apply(eEuDiff, c(1,2), dnorm, mean = 0, sd = 1, log = TRUE)
  #calculate log likelihood
  ell_vec=apply(ll_mat, 1, sum)


  print("nll")
  print(nll)

  nll<- -sum(ell_vec, na.rm=TRUE)
FicusBenji
  • 111
  • 5

0 Answers0