# [R] How might one write this better?

Tony Plate tplate at blackmesacapital.com
Tue Oct 5 16:59:30 CEST 2004

```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 = cont + ret * cont + ret * ret * cont

(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
   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))
 681.585
> # compare the two
> sum(cumprod(rev(c(1+ret[-1],1))) * rev(cont)) - asset[length(ret)+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