I was reading this seminal paper by Infanger. On page 40, Figure 11. was quite interesting. In particular I was interested in the top one, 19 Years and I wanted to reproduce this plot. To give some background: It's about utility maximization which should be solved by DP approach, i.e.
$$ \max_{x_t,0\le t\le T} E[u(W_T)]$$
where $u$ is a utility function and $W_T$ is the wealth at time $T$. We want therefore to maximize terminal wealth. For the picture he uses the following "quadratic downside risk" function
$$ u(W) = W - \frac{\lambda}{2}\max{(0,W_d-W)^2}$$
where $W_d$ is a target amound and $\lambda$ a scaling parameter. As wealth evolves via $W_{t+1} = W_t \cdot\langle x_t, R\rangle$ where $x_t$ are the allocation and $R$ the return, he writes down the Bellman equation of this problem:
$$V_{t}(W_t) = \max_{x_t}E[V_{t+1}(W_t\cdot \langle x_t, R\rangle)|W_t]$$ since I'm just interested in the final step, we have $V_T = u$ and the maximization problem I want to solve is $$V(W) = \max_{x}E[u(W\cdot \langle x, R\rangle)|W] $$ dropping the time $t$ index. I assumed (as Infanger did if I get him right), that $R$ are normally distributed. I use then Gauss-Hermite quadrature to approximate the expectation (see page 71 in this paper). $$V(W) = \frac{1}{\sqrt{\pi}}\max_{x}\sum_{i=1}^m w_i u(W(1+\hat{\mu}(x)+\sqrt{2}\hat{\sigma}(x)\cdot q_i)) $$ where $\hat{\mu} = \langle \mu, x\rangle$ and $\hat{\sigma} = \sqrt{\langle x, \Sigma x\rangle}$ and $w_i$ are the Gauss-Hermite weights and $q_i$ the corresponding nodes.
I've coded a very simple and not optimized version to see if I get the desired picture.
first a picture of the utility function:
utility <- function(w){
K <- 100000
temp <- K-w
temp[temp<=0] <- 0
return(w - 1000*temp^2)
}
x <- seq(90000,120000,1000)
y <- utility(x)
plot(x,y,type="l")
Now I just generated a sequences of $W$ and solved the above problem. The following code junk defines the covariance and expected return vector. The data is from the Infanger paper above.
wealth <- seq(50000,150000,5000)
mu <- c(0.108, 0.1037, 0.0949, 0.079, 0.0561)
cor <- matrix(c(1, 0.601, 0.247, 0.062, 0.094,
0.601, 1.0, 0.125, 0.027, 0.006,
0.247, 0.125, 1.0, 0.883, 0.194,
0.062, 0.027, 0.883, 1.0, 0.27,
0.094, 0.006, 0.194, 0.27, 1.0),
ncol=5, nrow=5,byrow=T)
std <- c(0.1572, 0.1675, 0.0657, 0.0489, 0.007)
temp <- std%*%t(std)
cov <- temp*cor
With this data at hand and a sequence of wealth (see above) I just run an optimization for each given wealth and store the solution (assuming no short selling). To solve the problem I used the Rsolnp package in R. It solves a minimization problem that's why I'm returning a $-1$ in the objective function below:
library(Rsolnp)
library(statmod)
obj <- function(x, currentWealth, mu, cov, r=0, nodes, weights){
drift <- sum((mu-r)*x)+r
cor <- sqrt(sum(x*(cov%*%x)))
term1 <- currentWealth*(1+drift)
term2 <- currentWealth*sqrt(2)*cor*nodes
return(-1/sqrt(pi)*sum(weights*utility(term1+term2)))
}
g_constraints <- function(x,currentWealth, mu, cov, r=0, nodes, weights){
return(sum(x))
}
x0 <- rep(0.25,length(mu))
weights <- gauss.quad(10,"hermite")$weights
nodes <- gauss.quad(10,"hermite")$nodes
solmat <- matrix(NA, ncol=length(mu),nrow=length(wealth))
for(i in 1:length(wealth)){
sol <- solnp(pars=x0, fun = obj,
eqfun = g_constraints,
eqB = 1,
LB = rep(0, length(mu)),
UB = rep(1, length(mu)),
currentWealth = wealth[i], mu = mu, cov = cov,
r = 0, nodes = nodes, weights = weights)
solmat[i,] <- sol$pars
x0 <- sol$pars
}
colnames(solmat) <- c("US Stock", "Int Stocks", "Corp Bonds", "Gvnt Bond", "Cash")
rownames(solmat) <- as.character(wealth)
However, I get a constant allocation where all money is invested in US Stocks. What's wrong with this and how do I get this chart from Infanger?
