[R] How might one write this better?

[Tony Plate]

The trick to vectorizing

 > asset <- numeric(T+1)
 > for (t in 1:T) asset[t+1] <- cont[t] + ret[t]*asset[t]

is to expand it algebraically into a sum of terms like:

asset[4] = cont[3] + ret[3] * cont[2] + ret[3] * ret[2] * cont[1]

(where the general case should be reasonably obvious, but is more work to 
write down)

Then recognize that this a sum of the elementwise product of a pair of 
vectors, one of which can be constructed with careful use of rev() and 
cumprod():

 > set.seed(1)
 > ret <- (rnorm(5)+1)/10
 > cont <- seq(along=ret)+100
 > asset <- numeric(length(ret)+1)
 > # loop way of computing assets -- final asset value is in the last 
element of asset[]
 > for (i in seq(along=ret)) asset[i+1] <- cont[i] + (1+ret[i]) * asset[i]
 > asset
[1]   0.0000 101.0000 214.9548 321.4880 508.9232 681.5849
 > # vectorized way of computing final asset value
 > sum(cumprod(rev(c(1+ret[-1],1))) * rev(cont))
[1] 681.585
 > # compare the two
 > sum(cumprod(rev(c(1+ret[-1],1))) * rev(cont)) - asset[length(ret)+1]
[1] 0
 >


At Sunday 05:35 AM 10/3/2004, you wrote:
>I am trying to simulate the trajectory of the pension assets of one
>person. In C-like syntax, it looks like this:
>
>daily.wage.growth = 1.001                 # deterministic
>contribution.rate = 0.08                  # deterministic 8%
>Wage = 10                                 # initial
>Asset = 0                                 # initial
>for (10,000 days) {
>     Asset += contribution.rate * Wage           # accreting contributions
>     Wage *= daily.wage.growth * Wage            # wage growth
>     Asset *= draw from a normal distribution    # Asset returns
>}
>cat("Terminal asset = ", Asset, "\n")
>
>How can one do this well in R? What I tried so far is to notice that
>the wage trajectory is deterministic, it does not change from one run
>to the next, and it can be done in one line. The asset returns
>trajectory can be obtained using a single call to rnorm(). Both these
>can be done nicely using R functions (if you're curious, I can give
>you my code). Using these, I efficiently get a vector of contributions
>c[] and a vector of returns r[]. But that still leaves the loop:
>
>   Asset <- 0
>   for (t in 1:T) {
>     Asset <- c[t] + r[t]*Asset
>   }
>
>How might one do this better?
>
>I find that using this code, it takes roughly 0.3 seconds per
>computation of Asset (on my dinky 500 MHz Celeron). I need to do
>50,000 of these every now and then, and it's a pain to have to wait 3
>hours. It'll be great if there is some neat R way to rewrite the
>little loop above.
>
>--
>Ajay Shah                                                   Consultant
>ajayshah at mayin.org                      Department of Economic Affairs
>http://www.mayin.org/ajayshah           Ministry of Finance, New Delhi
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://stat.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html