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