From: Ray Brownrigg <Ray.Brownrigg_at_mcs.vuw.ac.nz>

Date: Thu, 20 Mar 2008 13:58:27 +1300

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Thu 20 Mar 2008 - 01:06:14 GMT

Date: Thu, 20 Mar 2008 13:58:27 +1300

On Thu, 20 Mar 2008, Alberto Monteiro wrote:

> Jonas Malmros wrote:

*> > But what if I had a matrix, where the last column was filled with
**> > values first (again, a for loop), and the rest was filled by using a
**> > double loop?
**>
**> Killing one of those loops is quite simple. The other may be
**> harder.
**>
**> > OVal <- matrix(0, n+1, n+1)
**> >
**> > for(i in 0:n){
**> > OVal[i+1, n+1] <- max(Val[i+1, n+1]-K, 0)
**> > }
**>
**> This is straightforward: notice that most R arithmetic
**> functions that operate in scalars also operate in vectors,
**> matrices and arrays. For example, sin(pi) is zero, but
**> sin(c(0,pi/2,pi,3*pi/2)) is c(0, 1, 0, -1) (with some
**> rounding errors).
**>
**> The loop in _i_ is just to take the max of a column? So:
**>
**> OVal[1:(n+1), n+1] <- # this is a vector
**> max(Val[1:(n+1), n+1] - # this is another vector
**> - K, 0) # vector - scalar = vector, max(vector) = vector
**>
**> or, in one line:
**>
**> OVal[1:(n+1), n+1] <- max(Val[1:(n+1), n+1] - K, 0)
**>
**> You can drop the indices, as you are taking all of them:
**>
**> OVal[,n+1] <- max(Val[,n+1] - K, 0)
**>
**> > for(i in seq(n,1, by=-1)){
**> > for(j in 0:(i-1)){
**> > OVal[j+1, i] <- a*((1-p)*OVal[j+1, i+1]+p*OVal[j+2, i+1])
**> > }
**> > }
**>
**> The inner loop is simple, but somehow tricky.
**>
**> You are computing each j-th term as the same j-th term combined
**> with the (j+1)-th term. So, you take a combination of js in
**> the 1:i range and combine with js in the 2:(i+1) range... So:
**>
**> OVal[1:i, i] <- a*((1-p)*OVal[1:i, i+1] + p*OVal[2:(i+1), i+1])
**>
**> The outer loop (in i) probably can't be optimized.
**>
**> Alberto Monteiro
*

Unfortunately, neither of these loop removals gives the same answers as the original code :-( [Hence the request for reproducible code]. Try: TreeMy(K=0.5)

What does work though is the following:

TreeMy <-

function(S=50,K=50,sigma=0.4,r=0.1,Time=1,n=3){
dtime <- Time/n

u <- exp(sigma*sqrt(dtime)) d <- 1/u p <- (exp(r*dtime)-d)/(u-d) a <- exp(-r*dtime)

OVal <- matrix(0, n+1, n+1)

Val <- outer(0:n, 0:n, function(i, j) ifelse(i <= j, u^i*d^(j-i), 0))

OVal[, n+1] <- pmax(Val[, n+1]-K, 0)

for (i in n:1) {

OVal[1:i, i] <- colSums(matrix(a*c((1 - p), p)* OVal[, i + 1][rep(1:2, i) + rep(0:(i - 1), each=2)], ncol=i)) }

list("Stock Price"=Val, "Option Price"=OVal) }

Here the first loop is removed with a simple pmax() (the 'vector' form of max()). The second loop uses the trick of expanding the column operations into a vector using recycling then restructuring them back into a matrix. This provides a significant speed improvement.

**HTH
**

Ray Brownrigg

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Thu 20 Mar 2008 - 01:06:14 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Thu 20 Mar 2008 - 01:30:27 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*