Estaba leyendo este seminal papel por Infanger. En la página 40, la figura 11. era bastante interesante. En particular, me interesó la de arriba, 19 años, y quise reproducir este gráfico. Para dar algunos antecedentes: Se trata de la maximización de la utilidad que debe ser resuelto por el enfoque DP, es decir
$$ \max_{x_t,0\le t\le T} E[u(W_T)]$$
donde $u$ es una función de utilidad y $W_T$ es la riqueza en el momento $T$ . Queremos, por tanto, maximizar la riqueza terminal. Para la imagen utiliza la siguiente función de "riesgo cuadrático a la baja"
$$ u(W) = W - \frac{\lambda}{2}\max{(0,W_d-W)^2}$$
donde $W_d$ es un objetivo en torno a y $\lambda$ un parámetro de escala. A medida que la riqueza evoluciona a través de $W_{t+1} = W_t \cdot\langle x_t, R\rangle$ donde $x_t$ son la asignación y $R$ el retorno, escribe la ecuación de Bellman de este problema:
$$V_{t}(W_t) = \max_{x_t}E[V_{t+1}(W_t\cdot \langle x_t, R\rangle)|W_t]$$ ya que sólo me interesa el último paso, tenemos $V_T = u$ y el problema de maximización que quiero resolver es $$V(W) = \max_{x}E[u(W\cdot \langle x, R\rangle)|W] $$ dejar de lado el tiempo $t$ índice. Supuse (al igual que Infanger si lo entiendo bien), que $R$ se distribuyen normalmente. Utilizo entonces la cuadratura de Gauss-Hermite para aproximar la expectativa (véase la página 71 de este papel ). $$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)) $$ donde $\hat{\mu} = \langle \mu, x\rangle$ y $\hat{\sigma} = \sqrt{\langle x, \Sigma x\rangle}$ y $w_i$ son los pesos de Gauss-Hermite y $q_i$ los nodos correspondientes.
He codificado una versión muy simple y no optimizada para ver si consigo la imagen deseada.
primero una imagen de la función de utilidad:
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")
Ahora acabo de generar una secuencia de $W$ y solucionó el problema anterior. El siguiente código basura define la covarianza y el vector de retorno esperado. Los datos proceden del documento de Infanger mencionado anteriormente.
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
Con estos datos a mano y una secuencia de riquezas (ver arriba) simplemente ejecuto una optimización para cada riqueza dada y almaceno la solución (asumiendo que no hay ventas en corto). Para resolver el problema utilicé el programa Rsolnp
en R. Resuelve un problema de minimización por lo que devuelve un $-1$ en la función objetivo siguiente:
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)
Sin embargo, tengo una asignación constante en la que todo el dinero se invierte en acciones estadounidenses. ¿Qué pasa con esto y cómo puedo obtener este gráfico de Infanger?