title: Portfolio simulation using Cholesky decomposition
body: <p>I do a MC simulation of a portfolio with 4 assets (Bond, equity, equity, cash market)
  I use monthly steps and my simulation horizon is 10 years i.e. 120 steps. My final goal is to compute the yearly expected shortfall, i.e. taking the worst 5% of the Portfolio Returns.</p>
  
  <p>The simulation seems to be ok - at the first glance. However, I have the impression the drift dominates the process over time, so my expected shortfall is even positive for the long end. Also the expected shortfall eventually decreases when I increase the weights for equity. This is also true when I set the expected return for each asset to zero and hence the increased risk should drag the expected shortfall down. 
  I expect a bug in my code but can't see it. Any advice highly appreciated!</p>
  
  <pre><code>#maturity in years
  maturity &lt;- 10
  #Using monthly steps
  nsteps &lt;- maturity*12
  dt &lt;-  maturity / nsteps
  #number of assets
  nAssets = 4
  #number of simulations
  nTrails = 10000
  
  #expected return p.a. for each asset, stored in vector BM.mu
  BM.mu &lt;- rep(NA,nAssets)
  BM.mu[1] &lt;- 0.0072
  BM.mu[2] &lt;- 0.0365
  BM.mu[3] &lt;- 0.04702
  BM.mu[4] &lt;- 0.0005
  
  #defining variable size
  simulated.Returns   &lt;- array(NA,    dim = c(nsteps+1, nTrails, nAssets))
  cumulative.PortReturns  &lt;- matrix(rep(NA,nsteps*nTrails), nrow = nsteps, ncol = nTrails)
  ES  &lt;- rep(NA, maturity)
  
  #defining my monthly correlation and covariance matrix
  corr_matrix &lt;- matrix(c(1.000000000, -0.05081574, -0.07697585,  0.0051,
              -0.050815743,  1.00000000,  0.80997805, -0.3540,
              -0.076975850,  0.80997805,  1.00000000, -0.3130,
              0.005095699, -0.35365332, -0.31278506,  1.0000), nrow = 4, ncol = 4)
  
  cov_matrix &lt;- matrix(c(1.44e-04, -2.20e-05, -3.86e-05,  8.44e-08,
              -2.20e-05,  1.30e-03,  1.22e-03, -1.76e-05,
              -3.86e-05,  1.22e-03,  1.75e-03, -1.81e-05,
              8.44e-08, -1.76e-05, -1.81e-05,  1.90e-06), nrow = 4, ncol = 4)
  
  #defining my portfolio weights
  port.weights &lt;- c(0.72, 0.07, 0.07, 0.14)
  
  #performing cholesky decomposition
  R &lt;- (chol(corr_matrix))
  
  #generating standard-normal, random variables
  x &lt;- array(rnorm(nsteps*nTrails*nAssets), c(nsteps*nTrails,nAssets))
  
  #generating correlated standard-normal, random variables
  ep &lt;- x %*% R
  
  #defining the drift 
  drift &lt;- BM.mu - 0.5 * diag(cov_matrix) 
  
  #generating asset paths
  temp = array(exp(as.vector(drift %*% t(dt)) + t(ep *sqrt(diag(cov_matrix)))), c(nAssets,nsteps,nTrails))
  
  for(i in 2:nsteps) temp[,i,] = temp[,i,] * temp[,(i-1),] 
  
  #changing dimension of the array temp from dim(nAssets, nsteps, nTrails) to dim(nsteps, nAssets, nTrails)
  simulated.Returns &lt;- aperm(temp,c(2,1,3))
  
  #computing portfolio returns for each simulation (nTrails). To do this, each step is weighted with "port.weights"
  #Since I generate continuous returns, I first transform them into discrete, multiply with weights and then transform back into continuous.
  for (z in 1:nTrails) {
  for (i in 1:nsteps) cumulative.PortReturns[i,z] = log(1+((exp(simulated.Returns[i,,z]-1)-1) %*% port.weights))
  }
  
  #Finally I compute the monthly expected shortfall (5%-level) by taking the average of the 5% worst portfolio yields
  #I do steps of 12 as I calculate the ES at the end of each year
  z = 0
  for (i in seq(12, nsteps, by = 12 )) {
  z = z + 1 
  ES[z]   &lt;- mean(sort(cumulative.PortReturns[i,]) [1:(0.05*nTrails)])
  }
  
  #plotting a sample of simulated portfolio returns
  #library(QRM)
  plot(as.timeSeries(cumulative.PortReturns[,1:100]), plot.type = 'single')
  </code></pre>