expected.utility <- function( sigma, n.M, B, w, n.V, nominal.level, 
  c.star ) {

  p.c.star <- 1 - 2 * pnorm( - c.star / ( sigma * sqrt( 1 + 1 / n.M ) ) )

  S.threshold <- round( nominal.level * n.V )

  bonus.probability <- 1 - pbinom( S.threshold - 1, n.V, p.c.star )

  return( B * bonus.probability - 2 * c.star * w * n.V )

}

sigma <- 5

n.M <- 365

B <- 5000

w <- 0.5

n.V <- 365

nominal.level <- 0.95

c.star <- seq( 0, 20, length = 500 )

e.u <- expected.utility( sigma, n.M, B, w, n.V, nominal.level, c.star )

cbind( c.star, e.u )

plot( c.star, e.u, xlab = 'c*', ylab = 'Expected Utility',
  type = 'l' )

##################################################################

expected.utility.for.optim <- function( c.star ) {

  p.c.star <- 1 - 2 * pnorm( - c.star / ( sigma * sqrt( 1 + 1 / n.M ) ) )

  S.threshold <- round( nominal.level * n.V )

  bonus.probability <- 1 - pbinom( S.threshold - 1, n.V, p.c.star )

  return( - ( B * bonus.probability - 2 * c.star * w * n.V ) )

}

sigma <- 5

n.M <- 365

n.V <- 365

B <- 5000

w <- 0.5

nominal.level <- 0.95

print( c.star.optimal <- optim( 10, fn = expected.utility.for.optim, 
  method = "L-BFGS-B", lower = 5, upper = 15 )$par )

# 10.86412

1 - 2 * ( 1 - pnorm( c.star.optimal / 
  ( sigma * sqrt( 1 + 1 / n.M ) ) ) )

# 0.969982