Skip to content

How to remove double loop?

5 messages · Jonas Malmros, davidr at rhotrading.com, Alberto Monteiro +1 more

Original
Jonas Malmros
# 
Bill, Alberto, Gabor,
Thank you for answering my question. Now I learned about outer() function.
That was a straightforward example.

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?

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)
}

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])
    }
}

Even if I leave the first simple for loop as it is, is there a more
efficient way to program the double loop part now that OVal is used
within the function itself?
It is pretty easy to write for loops, but it is very hard to write
computationally optimal code. :-( Could you please help me with the
above one, if possible?
davidr at rhotrading.com
# 
Jonas,
Take a look at CRRBinomialTreeOption{fOptions} in which Diethelm Wuertz
uses a double loop,
but of course, he's optimized it to use only a vector instead of a
matrix!

David L. Reiner, PhD
Head Quant
Rho Trading Securities, LLC


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of Jonas Malmros
Sent: Wednesday, March 19, 2008 9:04 AM
To: r-help at r-project.org
Subject: [R] How to remove double loop?

Bill, Alberto, Gabor,
Thank you for answering my question. Now I learned about outer()
function.
That was a straightforward example.

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?

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)
}

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])
    }
}

Even if I leave the first simple for loop as it is, is there a more
efficient way to program the double loop part now that OVal is used
within the function itself?
It is pretty easy to write for loops, but it is very hard to write
computationally optimal code. :-( Could you please help me with the
above one, if possible?
Ray Brownrigg
# 
Read the posting guide (the bit about reproducible code):
+     OVal[i+1, n+1] <- max(Val[i+1, n+1]-K, 0)
+ }
Error: object "Val" not found
+     for(j in 0:(i-1)){
+         OVal[j+1, i] <- a*((1-p)*OVal[j+1, i+1]+p*OVal[j+2, i+1])
+     }
+ }
Error: object "p" not found

        
On Thu, 20 Mar 2008, Jonas Malmros wrote:

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Alberto Monteiro
    

    

      
      #
    

  


  

      
Jonas Malmros wrote:

            

      
      
      
    

        
Killing one of those loops is quite simple. The other may be
harder.

            

      
      
      
    

        
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)

            

      
      
      
    

        
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

  
  


  


      

    

    

      
      

        


  

        

      Ray Brownrigg
    

    

      
      #
    

  


  

      
On Thu, 20 Mar 2008, Alberto Monteiro wrote:

            

      
      
      
    

        
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