0

Apologies for the title. I could not think of how to title this...

(Also, I know this is probably a shite question, but hoping someone out there can help.)

I have the following mean vector and covariance matrix:

> mu0
       MSFT        AAPL 
0.001250251 0.001060690 
> sig0
             MSFT         AAPL
MSFT 1.275625e-04 3.334225e-05
AAPL 3.334225e-05 1.484212e-04

I have the following chunk of code (I know it is only 500 simulations. This is for brevity's sake. I will ramp this up when I figure out the code.):

wta = 0.562546911; wtm = 0.437453089; ct = -1000000*c(wtm,wta); lambda = .97; k = 20; 
M = 500; L1 = 0.0; l1 = rep(0,M); 
mu1 = mu0; sig1 = sig0; L1=0.0;
rando = rmvnorm(M, mu1, sig1)
for(i in 1:M){ 
  Xtd = t(rando[i,])
  mu1 = lambda*mu1+(1-lambda)*Xtd
  sig1 = lambda*sig0+(1-lambda)*t(Xtd)%*%Xtd; L1 = L1+t(ct)%*%t(Xtd)
  L1=L1+t(ct)%*%t(Xtd)
  for(j in 2:k){ 
    Xtd=rmvnorm(1,mu1, sig1)
    mu1=lambda*mu1+(1-lambda)*Xtd
    sig1=lambda*sig0+(1-lambda)*t(Xtd)%*%Xtd
    L1=L1+t(ct)%*%t(Xtd)
  } 
  l1[i]=L1
}

I hate this code. It is very slow, and perhaps more importantly, very ugly. I feel it is ripe for some sort of optimization using Reduce or do.call or something. But I cannot figure out how to write the functional for the Reduce function. To be perfectly honest, I do not fully understand Reduce. I have been trying to read through this, but struggling to apply it here.

The problem seems to be that I have several value which must 'accumulate' through the pulls of random numbers, so it is like my 'x' for Reduce is a list instead of a vector, which doesn't seem to make sense. I tried writing something like this:

runStuff <- function(thePrev, theNext, theLam, theCT, theSig){
  Xtd=rmvnorm(1, thePrev[1], thePrev[2])
  theNext[3]=thePrev[3]+t(theCT)%*%t(Xtd)
  theNext[1]=theLam*thePrev[1]+(1-theLam)*Xtd
  theNext[2]=theLam*theSig+(1-lambda)*t(Xtd)%*%Xtd
  theNext
}

And using it like this (which I fully expected to not work), and it does not work:

Reduce(x=list(mu1, sig1, L1), f=runStuff, 
       theLam = lambda, theCT = ct, theSig = sig0, accumulate=TRUE)

Any help with how I could use a functional in this scenario?

lukehawk
  • 1,423
  • 3
  • 22
  • 48
  • Can you replace your loops with an `sapply` statement? I've handled a similar Monte Carlo problem by doing as such. – Lloyd Christmas Nov 17 '16 at 15:05
  • The iterations of your outer loop seem to be fully independent.`Reduce` can only be useful for your inner loop, but it would be sufficient to [byte code compile](https://stat.ethz.ch/R-manual/R-devel/library/compiler/html/compile.html) your code if you find the `Reduce` syntax not to your taste. I see no benefit in the use of `lapply`/`sapply` here. What you should do to improve performance: Pre-allocate `l1 <- numeric(M)`. Then take everything out of the loop that doesn't change between iterations, such as `lambda*mu1` or `t(ct)`. Then study `help("crossprod")`. – Roland Nov 17 '16 at 15:17
  • And finally, consider if there might be a more efficient algorithm. – Roland Nov 17 '16 at 15:17
  • Do `Xtd=rmvnorm(M, mu1, sig1)`, use it as the first parameter of `Reduce` and pass the Xtd to `runStuff`. – Roland Nov 17 '16 at 15:20
  • I agree with Roland that sapply would not (in my understanding) do much here. I also agree that the outer loop is not the real issue. I should have been clearer, but I am really only concerned with the inner loop. Also, I already pre-allocated l1 with `l1= rep(0,M)`. – lukehawk Nov 17 '16 at 15:29
  • For runStuff, I cannot call `Xtd=rmvnorm(M, mu1, sig1)` since mu1 and sig1 change every iteration. I have to re-call rmvnorm with new mu1 and sig1 based on previous iterations. This is what I was trying to accomplish with `accumulate=TRUE`. – lukehawk Nov 17 '16 at 15:31
  • Then you shouldn't use `Reduce`. It's not more efficient than a byte-compiled `for` loop (because that's what it is internally) and won't be more readable in this case. – Roland Nov 18 '16 at 07:09

0 Answers0