utility function

util.r

####################################
# R code for a CRRA Utility function 
# Florian Oswald
# u(c) is defined over consumption c \in [-\infty,\infty]
# has automatic switch to log if CRRA = 1
# can produce diagnostic plot
# INPUT: x = vector of consumption values
# params   = list with 3 entries: CRRA(numeric), cutoff(numeric), diag.plot(bool)
#		 	 e.g. utility.params <- list('CRRA' = CRRA, 'cutoff' = cutoff, diag.plot=FALSE)
#			 'cutoff' is a number chosen by the user below which u(c) is quadratically approximated
# OUTPUT: list with approximated utility and gradient, floored utility, gradient and hessian and plot, if set
####################################

utility <- function(x,params){
	
	if (params$CRRA != 1){
		# CRRA utility (of floored consumption)
		floor.x                    <- x             # copy actual consumption (can be negative) into floored consumption
		floor.x[x < params$cutoff] <- params$cutoff # replace consumption below cutoff with cutoff
		diff.cons                  <- x - floor.x   # difference between actual and floored consumption

		tmp.x     <- floor.x^(1-params$CRRA)          # cons^(1-CRRA)
		du.dx     <- tmp.x/floor.x                    # u'(cons)
		ddu.dxx   <- -params$CRRA*du.dx/floor.x       # u"(cons)
		dddu.dxxx <- -(1+params$CRRA)*ddu.dxx/floor.x # u'''(cons)
		util      <- (tmp.x - 1)/(1-params$CRRA)      # u(cons)
	} else {
		# if CRRA==1, log utility
		floor.x                    <- x
		floor.x[x < params$cutoff] <- params$cutoff # replace consumption below cutoff with cutoff
		diff.cons                  <- x - floor.x   # difference between actual and floored consumption

		util      <- log(floor.x) # u(cons)
		du.dx     <- 1/floor.x    # u'(cons)
		ddu.dxx   <- -1/floor.x^2 # u"(cons)
		dddu.dxxx <- 2/floor.x^3  # u'''(cons)
	}

	# adjusting negative consumption values
	# quadratic approximation to utility and gradient, linear to hessian

	return.util <- util + du.dx*diff.cons + 0.5*ddu.dxx*diff.cons^2
	return.grad <- du.dx + ddu.dxx*diff.cons + 0.5*dddu.dxxx*diff.cons^2
	return.hess <- ddu.dxx + dddu.dxxx*diff.cons

	# make plot if diag.plot was set to TRUE
	if (params$diag.plot){
		par(mfcol=c(3,2))
		plot(x=x,return.util,type="l",main="approx. utility",xlab="consumption",ylab="utility")
		abline(v = params$cutoff,col="red")
		grid()
		plot(x=x,return.grad,type="l",main="approx. gradient",xlab="consumption",ylab="gradient")
		abline(v = params$cutoff,col="red")
		grid()
		plot(x=x,return.hess,type="l",main="approx. hessian",xlab="consumption",ylab="hessian")
		abline(v = params$cutoff,col="red")
		grid()
		plot(x=x,util,type="l",main="capped utility",xlab="consumption",ylab="capped utility")
		grid()
		plot(x=x,du.dx,type="l",main="capped gradient",xlab="consumption",ylab="capped gradient")
		grid()
		plot(x=x,ddu.dxx,type="l",main="capped hessian",xlab="consumption",ylab="capped hessian")
		grid()
		par(mfcol=c(1,1))
	}
	# return list of values
	return(list(utility=return.util,gradient=return.grad,hessian=return.hess,floor.util=util,floor.grad=du.dx,floor.hessian=ddu.dxx))
}