inst/BookEx/C14R1.R

In FRAPO: Financial Risk Modelling and Portfolio Optimisation with R

## Utility function
U1 <- function(x, mu, risk, lambda){
  lambda * x * mu - (1 - lambda) * risk * x^2
}
## Unscaled density
U1DU <- function(x, mu = 5, risk = 16, lambda = 0.5, nu = 1){
  exp(nu * U1(x = x, mu = mu, risk = risk, lambda = lambda))
}
## Constant of density
Z <- integrate(U1DU, lower = 0, upper = 1)$value
## Classical MC integration for expected value of scaled U1DU
MCsize <- 1e3
u1 <- runif(MCsize)
omega1 <- u1 * U1DU(u1) / Z * 100
omega1Hat <- mean(omega1)
## Graphical presentation of recursive means with bands
omega1RecMean <- cumsum(omega1) / (1:MCsize)
omega1RecErr <- sqrt(cumsum((omega1  - omega1RecMean)^2)) / (1:MCsize)
plot(omega1RecMean, type = "l",
     ylim = c(0, 100), ylab = expression(omega), xlab = "",
     main = "Recursive estimates for allocation to risky asset")
lines(omega1RecMean + 2 * omega1RecErr, col = "red", lty = 2)
lines(omega1RecMean - 2 * omega1RecErr, col = "red", lty = 2)
## Random draws from U1DU by accept-reject method
## Candidate distribution is a truncated Exponential with lambda  = 2 and k = 1.5
library(truncdist)
k <- 1.5
y2 <- rtrunc(MCsize, spec = "exp",
             a = 0.0, b = 1.0, rate = 2)
u2 <- runif(MCsize)
AcceptRejectIdx <- which(u2 <= (U1DU(y2) /
                                (k * dtrunc(y2, spec = "exp", a = 0, b = 1, rate = 2))))
omega2 <- y2[AcceptRejectIdx]
omega2Hat <- mean(omega2) * 100