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 [EMAIL PROTECTED] Department of Economic Affairs http://www.mayin.org/ajayshah Ministry of Finance, New Delhi
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
